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ABSTRACT  OF  THESIS 


FURTHER  DEVELOPMENT  AND  TESTING  OF  A  SECOND-ORDER  BULK 

BOUNDARY  LAYER  MODEL 


A  one-layer  bulk  boundary  layer  model  is  developed  following  earlier  work  by 
Randall  and  Moeng.  The  model  predicts  the  mixed  layer  values  of  the  potential 
temperature,  mixing  ratio,  and  u-  and  v-momentum.  The  model  also  predicts  the  depth  of 
the  boundary  layer  and  the  vertically  integrated  turbulence  kinetic  energy  (TKE).  The 
TKE  is  determined  using  a  second-order  closure  that  relates  the  rate  of  dissipation  to  the 
TKE.  The  fractional  area  covered  by  rising  motion  (a)  and  the  entrainment  rate  (E)  are 
diagnostically  determined. 

The  model  is  used  to  study  the  clear  convective  boundary  layer  (CBL)  using  data 
from  the  Wangara,  Australia  boundary  layer  experiment.  The  Wangara  data  is  also  used 
as  an  observation  base  to  validate  model  results.  A  further  study  is  accomplished  by 
simulating  the  planetary  boundary  layer  (PBL)  over  an  ocean  surface.  This  study  is 
designed  to  find  the  steady-state  solutions  of  the  prognostic  variables. 

The  model  clearly  illustrated  the  features  found  in  a  CBL.  The  diurnal  trend  of  the 
PBL  depth  was  accurately  reproduced.  This  included  rapid  growth  during  mid-morning, 
quasi-steady-state  conditions  during  the  afternoon,  and  an  evening  transition. 

In  the  ocean  study,  the  prognostic  variables  converged  to  their  equilibrium  values  at 
about  the  same  time.  This  is  in  contrast  to  an  earlier  study  using  similar  conditions  where 
the  adjustment  time  for  the  PBL  depth  was  considerably  longer  than  for  the  other 
prognostic  variables.  This  discrepancy  was  due  to  the  different  entrainment 
parameterizations  used  in  each  study.  In  the  ocean  study,  the  entrainment  rate  became 
very  large  during  the  initial  portion  of  the  simulation,  whereas  in  the  earlier  study  the 
entrainment  rate  remained  small  and  constant  throughout. 
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The  TKE  became  very  large  during  the  mid-morning  when  rapid  PBL  growth  was 
occurring.  This  large  TKE  indicated  that  the  PBL  was  very  turbulent  due  to  the  vigorous 
convection  that  was  taking  place.  The  fractional  area  covered  by  rising  motion,  a, 
reached  its  minimum  at  this  time;  a  further  indication  of  the  intense  convection. 

The  gradients  of  the  mean  potential  temperature  and  mean  mixing  ratio  were 
determined.  These  gradients  were  large  at  the  start  of  the  simulation  when  the  PBL  was 
unmixed.  The  gradients  decreased  rapidly  as  turbulence  mixed  the  PBL  during  mid- 
morning.  The  gradients  were  near  zero  in  the  afternoon  indicating  that  the  PBL  was  now 
well  mixed. 

A  two-layer  model  was  developed  to  address  the  problem  of  large  gradients  obtained 
in  the  one-layer  model.  This  model  produced  the  same  results  for  the  prognostic 
variables  as  the  one-layer  model.  The  gradients  determined  by  the  model  were  near  zero. 
The  mean  potential  temperatures  and  mixing  ratios  at  the  two  levels  in  the  model  were 
then  initially  perturbed  to  study  the  effects  of  varying  the  dissipation  time  scale.  A 
certain  range  of  values  of  the  model  parameter  related  to  the  dissipation  time  scale 
allowed  the  large  induced  gradients  to  approach  zero  in  a  reasonable  time. 

The  following  items  were  presented  for  the  first  time  in  this  thesis: 

(1)  A  positive  entrainment  rate  parameterization  which  assumes  a  balance  between 
buoyancy  production  and  dissipation  of  turbulence  kinetic  energy. 

(2)  A  negative  entrainment  rate  parameterization  that  allows  the  PBL  depth  to 
decrease  late  in  the  day  when  buoyancy  production  is  no  longer  sufficient  to  maintain  the 
turbulence. 

(3)  A  fully  implicit  finite  difference  equation  for  the  TKE  (when  the  entrainment  rate 
is  positive)  solved  as  a  cubic  equation.  The  square  of  the  solution  that  is  always  real  is 
assigned  to  the  TKE. 

(4)  Results  for  both  the  Wangara  and  Ocean  studies  showing  the  fractional  area 
cover<’d  by  rising  motion,  convective  mass  flux,  updraft  and  downdraft  properties  of  6 
and  q  at  the  surface  and  PBL  top,  dissipation  rates  of  9  and  q  at  the  surface  and  PBL  top. 
dissipation  time  scale,  and  gradients  of  6  and  q. 
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(5)  Results  and  comparison  for  the  Wangara  study  of  two  surface  bulk  transfer 
coefficients,  one  dependent  on  the  surface  velocity  and  the  other  on  the  turbulence  kinetic 
energy. 

(6)  A  two-layer  model  which  predicts  6  and  q  at  two  levels. 

(7)  Equations  that  determine  the  upward  turbulent  fluxes  of  0  and  q  in  the  interior  of 
the  PBL.  These  equations  are  used  to  obtain  0  and  q  in  the  two-layer  model. 

Richard  David  Krasner 
Department  of  Atmospheric  Science 
Colorado  State  University 
Fort  Collins,  CO  80523 
Summer  1993 


V 


ACKNOWLEDGMENTS 


I  would  like  to  express  my  sincere  appreciation  to  Dr.  David  Randall,  Dr.  Wayne 
Schubert,  Cindy  Carrick,  Douglas  Cripe,  Don  Dazlich,  Scott  Denning,  Jerry  Harrington, 
Ross  Heikes,  and  Debra  Youngblood  for  their  assistance  in  preparing  this  thesis. 

I  also  want  to  thank  my  wife.  Darken,  for  putting  up  with  me  during  the  past  18 
months.  Her  inspiration  kept  me  going  during  those  rough  moments. 

Support  for  this  research  was  provided  by  NASA  under  Grant  NAGl-1 137,  and  by  the 
Office  of  Naval  Research  under  Contract  NOO 14-9 1-J- 1422. 


VI 


TABLE  OF  CONTENTS 


1.  Introduction . 1 

La.  Statement  of  Problem . 1 

l.b.  Definition  of  Second-Order  Bulk  Boundary  Layer  Model . 2 

1. c.  Literature  Review . 3 

2.  Description  of  One-Layer  Model . 6 

2. a.  Equations . 6 

2.a.(l)  Conservation  of  Mass . 10 

2.a.(2)  Conservation  of  Momentum . 1 1 

2.a.(3)  Conservation  of  Potential  Temperature . 18 

2.a.(4)  Conservation  of  Moisture . 21 

2.a.(5)  Turbulence  Kinetic  Energy  (TKE)  Equation . 22 

2.a.(6)  Entrainment  Rate  Equations . 32 

2.a.(6)(a)  Positive  Entrainment . 35 

2,a.(6)(b)  Negative  Entrainment . 36 

2.b.  Initialization . 37 

2.C.  Top  Boundary  Conditions . 39 

2. d.  Surface  Boundary  Conditions . 40 

3.  One-Layer  Model  Time  Schemes . 41 

3. a  Surface  Heat-Moisture  and  Momentum  Flux  Parameterizations . 41 

3.b.  Conservation  of  Momentum . 43 

3.C.  Conservation  of  Potential  Temperature . 46 

3.d.  Conservation  of  Mixing  Ratio . 48 

3. e.  Turbulence  Kinetic  Energy . 49 

4.  Simulations . 51 

4. a.  Land  Simulation . 51 

4, b.  Ocean  Simulation . 51 

5.  One-Layer  Model  Prognostic  Results . 53 

5.  a.  Wangara  Experiments . 53 

5.a.(l)  Twenty-four  Hour  Simulation . 53 

5.a.(2)  Seventy-two  Hour  Simulation . 61 

5. b.  Ocean  Experiment . 62 

6.  One-Layer  Model  Diagnostic  Discussion . 70 

6. a.  Convective  Mass  Flux  Model . 70 

6.b.  Matching  Convective  Mass  Flux  with  Ventilation  and  Entrainment 

Mass  Flux . 72 

6.C.  Diagnostic  Equations  for  Me  and  a  Using  the  TKE . 73 

6.d.  PBL  Interior  Diagnostics . 75 

6.e.  Surface  Transfer  Coefficient  Using  TKE . 78 

6. f.  Richardson  Number  and  Limits . 78 

7.  One-Layer  Model  Diagnostic  Results . 80 

7.  a.  Wangara  Results  for  the  Fractional  Area  Covered  by  Rising  Motion . 80 

7.b.  Wangara  PBL  Interior  Results . 86 

l.c.  Wangara  Surface  Transfer  Coefficients . 91 

7.d.  Wangara  Calculation  of  Richardson  Number  and  Limits . 93 


vii 


7.e.  Ocean  Experiment  Fractional  Area  Covered  by  Rising  Motion 

Results .  ytj 

7. f.  Ocean  Experiment  PEL  Interior  Results . 102 

8.  Description  of  Two-Layer  Model . 107 

8.  a.  Two-Layer  Potential  Temperature  and  Mixing  Ratio  Equations . 108 

8. b.  Two-Layer  Model  Diagnostics . 113 

9.  Two- Layer  Model  Results . 114 

9.  a.  Two-Layer  Prognostic  Results . 114 

9.b.  Two-Layer  Diagnostic  Results . . . 117 

10.  Summary  and  Conclusions . 122 

References . 129 


viii 


1.  Introduction 


La.  Statement  of  Problem 

The  planetary  boundary  layer  (PBL)  exerts  a  significant  influence  on  the  earth’s 
weather  and  climate.  The  large-scale  atmosphere  feels  the  effects  of  the  PBL  during  the 
development  and  growth  of  thunderstorms  with  relatively  short  time  scales  («  1  hour), 
and  over  the  entire  time  scale  spectrum  including  long  periods  (lO’s  -  lOO’sof  years) 
when  global  climatic  change  occurs.  The  predominant  source  of  energy  to  drive  the 
general  circulation  is  the  ocean.  Surface  fluxes  of  heat,  moisture,  and  momentum  over 
the  oceans  are  transported  to  the  free  atmosphere  through  the  PBL.  These  fluxes  play  a 
vital  role  in  transforming  the  earth’s  climate  over  time.  The  turbulent  eddies  in  the  PBL 
are  the  means  by  which  this  energy  is  transmitted  from  the  ocean  to  the  free  atmosphere 
where  it  interacts  with  the  general  circulation.  Since  the  PBL  is  intimetely  tied  to  the 
evolution  of  the  climate,  an  accurate  representation  of  the  PBL  is  required  to  correctly 
make  predictions  of  the  future  climate  using  a  general  circulation  model.  The  only  way 
to  accomplish  this  is  to  develop  a  model  that  predicts  the  state  of  the  PBL  using 
parameterizations. 

The  surface  fluxes  of  heat,  moisture,  and  momentum  are  not  the  only  important 
parameters  to  consider.  The  fluxes  of  these  quantities  over  the  entire  depth  of  the  PBL 
should  also  be  included  because  they  affect  the  general  circulation.  The  total  water  vapor 
in  the  PBL  represents  the  latent  heat  available  to  drive  the  general  circulation.  Another 
parameter,  the  PBL  depth,  besides  denoting  the  amount  of  mass  contained  within  the 
PBL,  gives  insight  into  whether  clouds  are  present.  As  the  PBL  depth  increases,  moisture 
can  penetrate  higher  into  the  atmosphere,  eventually  reaching  the  lifting  condensation 
level  where  clouds  will  form.  Finally,  one  can  obtain  some  information  about  the 
fractional  cloud  amount  by  predicting  the  fractional  area  in  the  PBL  that  rising  motion 
covers.  Clouds  play  a  key  role  in  the  climate  because  they  affect  the  radiation  budget  by 
reflecting  solar  radiation  that  the  earth's  surface  would  otherwise  absorb,  therefore, 
knowledge  of  the  amount  of  cloud  coverage  is  crucial  to  climate  prediction. 
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l.b.  Definition  of  Second-Order  Bulk  Boundary  Layer  Model 

The  bulk  boundary  layer  model  presented  is  further  development  of  the  work 
completed  by  Randall,  Shao,  and  Moeng(1992).  The  model  is  1-dimensional  and 
employs  a  second-order  turbulence  closure,  and  a  "bulk"  approach  to  parametrically 
represent  boundary-layer  structure.  Prediction  equations  are  used  to  compute  the 
boundary  layer  depth  (Ap^.  in  terms  of  pressure),  mixed  layer  values  of  u-  and  v- 
momentum  (um  and  Vm),  potential  temperature  (Bm),  mixing  ratio  (qm),  and  turbulence 
kinetic  energy  (em).  This  model  also  diagnoses  the  entrainment  rate  (E)  and  fractional 
area  covered  by  rising  motion  (o),  which  allows  the  determination  of  fractional  cloud 
amounts.  Figure  l.b.l  depicts  the  domain  of  the  model.  Subscript  S-  denotes  the  earth's 
surface,  S  the  top  of  the  ventilation  or  surface  layer,  B  the  base  of  the  entrainment  layer 
or  the  top  of  the  PBL,  and  B+  the  top  of  the  entrainment  layer  or  the  level  just  above  the 
top  of  the  PBL.  Both  the  surface  and  entrainment  layers  are  infinitesimally  thick 
(indicated  by  the  stippling  in  figure  l.b.l). 


Figure  l.b.l:  Domain  of  Bulk  Boundary  Layer  Model. 
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l.c.  Literature  Review 


Deardorff  (1974a)  completed  a  three-dimensional  numerical  study  of  the  heated  PEL 
where  he  determined  the  mean  structure  and  height  of  the  PEL.  His  study  utilized  data 
from  Day  33  of  the  Wangara  Experiment:  Boundary  Layer  Data  (Clarke  et  al.,  1971). 

His  model  included  grid-volume  averaged  equations  for  the  momentum,  potential 
temperature,  and  mixing  ratio.  The  potential  temperature  equation  included  a  term  fur  the 
temperature  change  due  to  the  divergence  of  the  long-wave  radiative  flux.  Deardorff  also 

utilized  subgrid  transrort  equations  (-^m/m',  ^u'd',  ^d'q' 

at  ot  ot  at  at  at 

for  the  subgrid  Reynolds  stresses.  He  assumed  that  the  terrain  was  flat  and  the  surface 

temperature  and  surface  roughness  were  horizontally  homogeneous. 


The  surface  momentum  flux  was  prescribed  using  the  surface-layer  formulations  of 
Businger  et  al.  (1971)  and  the  surface  layer  integrals  of  Paulson  (1970).  The  surface  heat 
and  moisture  fluxes  were  computed  using  the  subgrid  transport  equations.  The  model 
was  initialized  using  data  beginning  at  0900L  on  Day  33. 

Deardorff  s  model  overestimated  the  calculated  rate  of  growth  of  the  mixed  layer 
height  between  1200-1500L.  He  attributed  this  overe'timation  to  the  lack  of  large-scale 
vertical  motion  in  the  model.  Thus,  the  model  mixed  layer  height  was  also  overestimated 
during  the  afternoon.  At  1200L  the  model  predicted  1030  meters,  while  the  actual  height 
was  1010  meters.  At  1500,  the  model  height  was  1400  meters  and  the  actual  height  was 
only  1200  meters.  The  maximum  model  height  was  1500  meters  at  1800L  and  the 
maximum  actual  height  vas  1280  meters. 


Yamada  and  Mellor  (1975)  completed  a  simulation  of  the  diurnally  varying 
planetary  bouna,.  ’ayer  and  compared  it  with  Days  33-34  of  the  Wangara  data.  Their 
model  differed  from  Deardorff  s  by  using  ensemble  mean  closure  instead  of  subgrid-scale 
closure.  They  used  their  level  3  model  which  required  the  solution  of  only  2  out  of  10 
differential  equations  for  turbulence  moments  (turbulence  kinetic  energy  and  temperature 
variance). 


Yamada's  model  underestimated  the  height  of  the  PEL  by  about  300  meters  during 
the  afternoon  hours.  This  discrepancy  was  due  to  the  uncertainty  in  the  observed  values 
of  the  mean  vertical  wind.  They  concluded  that  accurate  data  for  the  thermal  wind  and 
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mean  vertical  wind  are  necessary  to  obtain  realistic  simulations  for  the  mean  winds  and 
temperature. 

Suarez  et  al.  (1983)  developed  a  parameterization  for  the  PBL  to  be  used  in  the 
UCLA  General  Circulation  Model.  Their  work  provides  some  of  the  basis  for  this  thesis. 
The  parameterization  used  a  mixed-layer  approach  where  the  discontinuities  in 
temperature  and  moisture  at  the  top  of  the  PBL  were  modeled  using  jumps.  They  utilized 
a  modified  a-coordinate  and  bulk  equations  for  the  PBL. 

Based  on  the  earlier  work  by  Deardorff  (1970, 1972,  and  1974a,b),  the  PBL  depth 
was  determined  using  a  prognostic  equation.  Deardorff  showed  that  entrainment  was 
vitally  important  in  the  determination  of  the  PBL  depth  and  that  the  PBL  depth  affected 
the  bulk  Richardson  number,  which  in  turn  was  related  to  the  stability  dependence  of  the 
surface  transfer  coefficients. 

The  modified  G-coordinate  was  used  to  allow  the  varying  depth  of  the  PBL  to  be 
included  into  the  GCM.  In  the  conventional  G-coordinate  system,  surfaces  follow  the 
earth's  topography.  In  the  modified  G-coordinate  system,  the  earth's  surface  and  the  PBL 
top  are  both  coordinate  surfaces.  Th*s  allowed  the  PBL  to  be  more  effectively  coupled  to 
the  large-scale  dynamics.  Also,  the  detailed  structure  that  occurs  at  the  PBL  top  where 
the  jumps  are  did  not  have  to  be  resolved  by  the  GCM  grid  since  the  structure  was  at  the 
interface  of  the  two  lowest  GCM  layers  in  the  modified  G-coordinate  system. 

The  bulk  equations  at  the  PBL  top  related  the  flux  of  a  quantity  to  the  product  of  the 
entrainment  rate  and  the  change  (jump)  of  the  quantity  across  the  top  interface.  Surface 
bulk  formula  were  determined  using  similarity  theory  formulated  by  Businger  et  al. 
(1971),  just  as  Deardorff  used.  Since  the  fluxes  at  the  PBL  top  use  the  entrainment  rate, 
it  had  to  be  parameterized  in  terms  of  the  prognostic  variables.  This  parameterization 
was  based  on  separating  the  buoyancy  and  shear  terms  into  positive  and  negative 
production.  The  entrainment  rate  was  then  given  by  positive  production  minus  negative 
production  minus  dissipation. 

Randall  et  al.  (1992)  developed  a  .second-order  bulk  boundary-layer  model.  This 
model  matches  the  fluxes  associated  with  the  convective  mass  flux  with  the  surface  or 
ventilation  mass  flux  and  the  entrainment  mass  flux.  The  model  also  provides  the  first 
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physically  based  means  to  determine  the  fractional  area  covered  by  rising  motion,  o. 
Finally,  the  model  allows  the  "well-mixed"  assumption  to  be  relaxed. 


This  thesis  is  a  continuation  of  the  work  done  by  Randall  et  al.  (1992).  A  new 
entrainment  parameterization  is  introduced  based  on  the  sums  of  the  buoyancy  and  shear 
production  terms.  A  surface  transfer  coefficient  is  calculated  based  on  the  predicted 
turbulence  kinetic  energy.  Results  are  shown  using  the  diagnostics  developed  in 
Randall's  bulk  model  including  the  convective  mass  flux  and  the  fractional  area  covered 
by  rising  motion.  Finally,  a  two-layer  model  is  developed  and  tested. 
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2.  Description  of  One-Layer  Model 


2.a.  Equations 

This  chapter  provides  detailed  derivations  of  the  prediction  equations  used  by  the 
model.  Table  2.a.l  provides  a  summary  of  the  final  equations  used,  the  assumptions 
made  in  the  derivations,  and  the  boundary  conditions  applied  to  simplify  and  solve  the 
prognostic  equations. 

Tensor  notation  is  used  throughout  the  derivations  with  the  subscript  j  referring  to 
one  of  the  components  of  momentum  or  Cartesian  coordinates  (u=ui,  v=  U2,  w=U3,  x=xi, 
y=X2,  z=X3,  where  subscripts  i,  j,  or  k  are  used  with  the  values  1, 2,  or  3). 
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Table  2.a.l:  Model  Equation  Summary. 


Equation 


Assumptions  Top  B.C. 


Surface  B.C. 


Conservation  of  Mass 
(2.a.l) 


u-Momentum 

(2.a.2) 

_  E{us^-u„) 


p^M  v„  u 

§  <^C  fftOfft  ftt  i 

P»,diZ„ 

I  m  m 


Pm^m 


•  Boussinesq 

•  Shallow 
convection 

•  Viscous  effects 
ignored 

d  j  ^  f 

•  —  and  —  of 


JV'n  ’g) 

perturbation 
quantities  zero 

•  Mixed  layer 
values  of  V,  p, 
and  Az  used 

•  Horizontal 
homogeneity 

v-Momentum 

•  Same  as  u- 

(2.a.3) 

momentum 

equation 

1 

1 

+ 

1 

E 

•  Turbulence 
vanishes 

•  Momentum 
flux  vanishes 

•  Momentum 
flux  at  level  B 
vanishes  if  E<0 


•  Same  as  u- 

momentum 

equation 


•  Horizontal 
homogeneity 

•  No  mass 
crosses  earth's 
surface 

•  Louis  (1979) 
surface 

momentum  flux 
parameterization 


•  Same  as  u- 

momentum 

equation 


Table  2.a.l:  Model  Equation  Summary  (continued). 


Equation 

Assumptions 

Top  B.C. 

Surface  B.C. 

Potential  Temperature 
(2.a.4) 

Psfc^  l^m l^ea/— moijf.  ^  \ 

•  Same  as 
momentum 
equations 
except: 

•  Molecular 
conduction  and 
radiation 
divergence 
ignored 

•  Mixed  layer 
value  of  0  also 
used 

•  No  phase 
changes  of  water 

•  Turbulence 
vanishes 

•  Heat  flux 
vanishes 

•  Heat  flux  at 
level  B  vanishes 
ifE<0 

•  Same  as 
momentum 
equations 
except: 

•  Louis  (1979) 
surface  heat- 
moisture  flux 
parameterization 

0.74P.AZ. 

Moisture 

(2.a.5) 

•  Same  as 
momentum 
equations 
except: 

•  Source-sink 
term  assumed  to 
be  a  mean 
forcing  and 
equal  to  zero 

•  Molecular 
diffusion 
ignored 

•  Mixed  layer 
value  of  q  also 
used 

•  No  clouds 

•  Turbulence 
vanishes 

•  Moisture  flux 
vanishes 

•  Moisture  flux 
at  level  B 
vanishes  if  E<0 

•  Same  as 

potential 

temperature 

0.74p,Az, 
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Table  2.a.l:  Model  Equation  Summary  (continued!. 


Equation 

Assumptions 

Top  B.C. 

Surface  B.C. 

TKE(E>0) 

•  Same  as 

•  TKE  vanishes 

•  No  mass 

(2.a.6) 

momentum 

crosses  earth's 

^'"=  (S+B  Ee„ 

equations 

except: 

•  Vertical 
turbulent  flux  of 

surface 

TKE  at  top  equal 

•  Vertical 

•  Viscous 

to  value  at 

turbulent  flux  of 

D) 

dissipation  term 

surface 

TKE  at  surface 

approximated 

•  Pressure 

equal  to  value  at 
top 

•  Gravity  waves 

correlation  at  top 

neglected 

zero  (due  to 

•  Pressure 

neglect  of 

correlation  at 

•  Mixed  layer 

gravity  waves) 

surface  zero 

value  of  TKE 

(turbulence 

also  used 

•  Vector 
momentum. 

vanishes) 

•  constant 

heat,  and 
moisture  fluxes 

•  Vector 
momentum  flux 

in  surface  layer 

vanish 

constant  in 

and  is  parallel  to 
wind 

surface  layer  and 
is  parallel  to 
wind 

•  decreases 
linearly  in 

•  Momentum 

entrainment 

vanishes 

layer  and  is 
parallel  to  wind 

•  Horizontal 
homogeneity 

•  Louis  (1979) 

c.TFq 

1 

surface  heat- 
moisture  flux 
i  parameterization 

•  Ratio  of  Fsv 
and  pressure 
approximately 
linear 

2 

•  o  =  — 

ai 

(second  order 
closure) 
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Table  2.a.l:  Model  Equation  Summary  (continued). 


Equation 

Assumptions 

TKE  (E<0) 

(2.a.7) 

^  =  -^—weighti  Bo  + 

4'. 

S.-D) 

•  Same  as  TKE 
(E>0)  except: 

•  Production  of 
TKE  weighted 
based  on 
contribution  due 
to  local  change 
and  entrainment 
production 

Positive  Entrainment 
(2.a.8) 

•  No  clouds 

•  Turbulence 
required  to  mix 
newly  entrained 
air 

•  Balance 
between  buoyant 
production  and 
dissipation  of 
TKE 

Negative  Entrainment 
(2.a.9) 

^^(1- weigh,). 

S.-D) 

•  No  clouds 

•  Turbulence 
required  to  mix 
newly  entrained 
air 

•  Local  rate  of 
change  of  TKE 
neglected 

Top  B.C. 


Surface  B.C. 


Same  as  TKE 
(E>0)  except: 

•  Vector 
momentum, 
heat,  and 
moisture  fluxes 
vanish  at  level  B 
and  above  since 
E<0 


•  Same  as  TKE 
(E>0) 


2.a.(l)  Conservation  of  Mass 


The  mass  conservation  equation  is  expressed  by 


dp  ^  d(puj)  _  ^ 
dt  dXj 


(2.a.(l).l) 
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This  equation  is  not  used  directly  (the  form  of  this  equation  listed  in  Table  2.a.l  is  used  to 
predict  the  boundary  layer  depth),  but  is  used  to  derive  the  other  conservation  equations. 

2.a.(2)  Conservation  of  Momentum 


The  Navier-Stokes  equation  is 


du:  o  n  {^+ li)  d  dU: 

IT  =  IT-  (2.a.(2).  1) 


dx, 


p  BXi  p  dx  j  p  dXi  dxj 


where  f  is  the  coriolis  parameter  (2Qsin(j>),  is  the  alternating  unit  tensor  where 


^ij3  - 


+lfor  i  =  1  and  j  =  2 
-Ifor  i  =  2  and  j  =  1 
Ofori  =  j, 


p  is  the  dynamic  viscosity  coefficient,  p  is  the  density,  and  A,=-2p/3.  The  incompressible 
form  of  (2.a.(2).l)  is 


dU;  dU;  5.  ,  1  dp  U  d^U: 


(2.a.(2).2) 


Equation  (2.a.(2).2)  is  multiplied  by  the  density  (split  into  mean  and  perturbation  parts 
when  multiplied  by  g,  constant  otherwise;  the  Boussinesq  approximation)  to  give 


=  -pSu8  -  S,38  +  Pf^ijsUj  “  ^  ^  ^ 


(2.a.(2).3) 


Next,  (2.a.(2).3)  is  divided  by  the  mean  density  to  obtain 


dUi  du-  ™  d  ^  .  1  dp  d^u- 


(2.a.(2).4) 


where  v=  pi p  is  the  kinematic  viscosity.  For  shallow  convection  ^ 

P 


approximation  is  then  applied  to  (2.a.(2).4)  which  yields 


This 
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Now,  equation  (2.a.(2).5)  is  multiplied  by  p  and  is  added  to  Ui  multiplied  by  the 
continuity  equation  to  get 


dt 


—  BUi  ^  -  Bu,  B(  pUj  j  _ 

^  ^  ‘  Bt  ^  ^Bx^  ‘  Bx^  ^ 


g- 

^  > 

g 

_ .  -1  Bp  -  B^u. 


Equation  (2.a.(2).6)  is  then  put  into  flux  form. 


Bt  BXj 

^  1  Bp  _  B^U: 


*■  ft 


g 


V  / 


+  Pf^ip^j  - 


(2.a.(2).6) 


(2.a.(2).7) 


Only  the  horizontal  momentum  equations  (i=l,  2)  are  considered,  the  geostrophic  wind 
definition  is  used,  and  viscous  effects  are  ignored.  Then,  (2.a.(2).7)  becomes 


Bt  Bxj 


(2.a.(2).8) 


where  the  g  subscript  denotes  a  geostrophic  wind  component.  The  wind  components  are 
next  split  into  mean  and  perturbation  parts  to  give 


B{  pu-  +  pm'  )  _  B(  pUiUj  +  pUiUj  +  pi/iUj  +  pU^Uj )  ^ 

^  (2.a.(2).9) 

The  split  is  not  necessary  for  the  geostrophic  component  since  this  component  is  a 
constant.  Equation  (2.a.(2).9)  is  then  simplified  using  the  method  of  Reynolds  averaging. 


Consider  an  instantaneous  quantity.  A,  which  is  split  into  a  mean  and  perturbation 
component  (A  and  A').  The  mean  component  of  A  represents  either  the  time,  space,  or 
ensemble  average  of  A,  and  the  perturbation  component  represents  positive  or  negative 
deviations  from  this  average.  If  A  =  A  +  A',  then  (a)  =  |a  +  A'j,  or  A  =  A  +  A'.  The 

last  equality  can  only  be  true  if  A'  =  0  This  just  states  that  the  sum  of  positive  deviations 
from  the  mean  equals  the  absolute  value  of  the  sum  of  the  negative  deviations,  thus  the 
net  sum  of  the  deviations  is  zero.  Reynolds  averaging  is  accomplished  by  applying  the 
above  result  to  quantities  split  into  mean  and  perturbation  parts.  Stull  (1991)  provides  a 
detailed  discussion  of  Reynolds  averaging.  Equation  (2.a.(2).9)  then  becomes 


djpUj)  _  S{pu,Uj  + 
dt  dXj 


(2.a.(2).10) 


The  u-component  of  (2.a.(2).10)  is 


d[pu)  _  d[ pu u)  S{ pu  v)  d{ pu  w) 


dt 


Bx 


dy 


Bz 


^P^^)  . 

Bx  By 


Bz 


(2.a.(2).ll) 


Horizontal  derivatives  of  perturbation  quantities  are  neglected  because 

.  .  *  u  A(horizontal flux)  A(vertical flux)  ^ 

Ax  «  Ay  »  Az,  thus - - — - - -  « - - - .  Equation  (2.a.(2).l  1), 


Ax,  Ay 


Az 


without  these  terms,  is  next  vertically  integrated  from  the  lower  surface  in  the  boundary 
layer,  zs-,  to  Just  above  the  boundary  layer  top,  zb+  which  yields 


Zit  -XT7 


'TB{pi^ 


f  f 

J  Bt  J 


z=rs- 


Z=Zs_ 


Z=Js.  -\f  — 


Bx 


dz-  I 


B{  pu  v) 


dz- 


B{  puw)  ^  ^  pi^  w  ) 


I  I 


2=JS- 


Bz 

J  P/(v-vJdz. 


!=^S- 


Bz 


dz  + 


»=2s- 


(2.a.(2).12) 


Equation  (2.a.(2).12)  is  then  transformed  using  Leibniz’s  rule  to  give 
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^  j  pudz=p{t.  Zs^)u{t.  ^s-)^ 

^~^S~ 

-—  J  puudz+p{t,  Zg^)u{t,  Zgju(t,  ZgJ-^ 


J=ZS- 


^  ^  r 

p(^  v)«{^  Zs_)M(r.  v)-^--^  J  pMvyz+ 


p(^  2b+)^- 


p(r,  Zj_)u(?,  Zs_)v(t.  Zs_) 


^s- 

dy 


P(^  ^bM^’  ^B*)^t’  ^fl+)  + 

p(^  Z5_)u(r.  Zj_)iv(/,  Z5_)-p(f.  ZgJi/M/(f,  z^^)  + 
p{/.  Zs_)i/W[t,  Zs_)  +  pJ(v„-v^)Az^, 


(2.a.(2).13) 


where  the  mixed  layer  values  (denoted  by  subscript  m)  of  v  (the  v-component  of  the 
wind),  the  density,  and  the  boundary  layer  depth  have  been  used  to  simplify  the 
integration.  Since  turbulence  vanishes  above  the  boundary  layer,  p(t,  ZB+)u'w'(t,  z^^.)  is 

zero. 


Next,  the  terms  at  height  B+  and  at  height  S-  are  combined  to  give 


—  J  pudz  =  it 


dt 


\  dt  dx  dy 


yj 


u 


—{ dz^  ^dz^  ^dz^ 
p  —^+u-r^+v—:r — w 


^  dt  dx  dy 

-J  2=2^  ^  25= 

d  (  d  t 


/  J 


(2.a.(2).14) 


-—  J  puudz-—  j  puvdz  + 
(  P^)s_  +  PJ(  -  \ 


where  the  function  notation  has  been  dropped.  The  terms  inside  the  square  brackets 
represent  the  entrainment  mass  flux,  E,  across  the  B+  and  S-  surfaces  respectively.  The 
S-  surface  is  the  earth’s  surface,  where  the  entrainment  mass  flux  is  zero  (the  individual 
terms  are  not  necessarily  zero,  but  their  sum  must  be  zero  since  no  mass  can  cross  the 
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earth’s  surface).  Equation  (2.a.(2).14)  then  simplifies  to 


—  J  pudz=  -V  •  J  upvdz  +  (pi/W)^  +  pj(v„  -  v^)Az„,  (2.a.(2).15) 

2=ZS-  *=Zs- 


where  the  second  term  on  the  right  hand  side  is  a  combination  of  the  fourth  and  fifth 
terms  in  equation  (2.a.(2).14).  Now,  the  first  term  and  the  divergence  term  are  integrated 
using  mixed  layer  values  as  was  done  for  the  coriolis  term  to  give 


or 


-V  •  [ +  (pi/ m/ )^_  + 
Pm/(v„-vjdz„, 


Pm4»2„  ^  -  U„  )  + 

{pu^\_  +  PJ{v„-v^)Az^, 


(2.a.(2).16) 


The  derivative  in  the  second  term  on  the  right  hand  side  of  the  lower  equation  of 
(2.a.(2),16)  can  be  expressed  in  terms  of  the  entrainment  mass  flux  into  the  top  of  the 
boundary  layer.  This  is  shown  by  Figure  2.a.(2).l. 
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(1)  Entrainment  brings  mass 
into  PEL  top  at  local  point 

(Pm 


(2)  Mass  flux  converges 
at  local  point 


(Pm^)l 

Area  of  Area  of 

more  mass  less  mass 

(3)  Area  containing  more  mass  is 
advected  into  local  point  by  the 
velocity 

Figure  2.a.(2).l:  Processes  Which  Cause  Local  Change  in  Mass. 

The  local  mass  flux  changes  due  to;  (1)  the  entrainment  of  mass  into  the 

PEL  top  (E),  (2)  horizontal  convergence  of  mass  flux  (-p„A2„(V  •  and  (3) 

horizontal  advection  of  mass  flux  (-v„  Therefore, 

=  E- P.Az„(V •  -  v„  •  V(p„Az„),^_^,.  Substituting  this  into  the  last 

equation  of  (2.a.(2).16)  gives. 

P«^z„  ^  =  EUs^-u„[e-  p„Az„ (V •  v„ )  -  v„  •  V(p„ Az„ )] 

-V  •  (m„P„v„ Az„)  +  (p;iV)^_  +  p„  f (v„  -  vJAz„ , 

or  using  the  vector  identity  u„V  •  {p„v„/k„)  =  M„p„Az„(V  •  v„)  +  •  V(p„Az„), 

Pm^^n,  ^  -  Eu„  +  M„ V  •  (p„ Az„v^ ) 

-V  •(m„P„v„Az„)  +  (pmV')^_  +  pj(v„  -  v,)Az„, 

or 
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V  •  (p„v„/lz„ )  -  M„  V  •  {p„v„Az„, )  - 

•  Vm„  +  +  pj(v„  -  )Az„  , 

^  _ 

-^  =  E{ub,-u„)- P^Az^v^  •  V«„  +  (pu'w\_ 

+P./(v„-v,)4z„. 


Horizontal  homogeneity  is  assumed  in  the  last  equation  of  (2.a.(2).17).  This 
eliminates  the  second  term  on  the  right  hand  side.  Then  this  equation  is  divided  by  the 
mixed  layer  density  and  boundary  layer  depth  to  give 


(2.a.(2).18) 


The  local  rale  of  change  of  the  mixed  layer  u-momentum  component  is  due  to  the 
entrainment  flux  of  momentum  through  the  boundary  layer  top,  the  surface  flux  of 
momentum,  coriolis-pressure  gradient  effects,  and  momentum  divergence.  The  v- 
component  equation  is  derived  in  the  same  manner  and  is 


^  Pn,^n,  Pm^n,  ^  ' 


(2.a.(2).19) 


Chapter  3  provides  details  of  the  parameterization  of  the  surface  momentum  fluxes. 


The  equations  with  the  included  parameiv  -ization  are 


^  E{ub,-u„)  p^.a^F„.„.KK  /  _  V 


and 


^'m  _  ^{^’b+  _  Psfc^  ^mom  _  rL  _  , ,  \ 

St  P^Az^  '  PA 


(2.a.(2).20) 


(2.a.(2).21) 
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where  is  the  surface  air  density,  a^  is  a  drag  coefficient,  is  an  empirical 
function  which  is  dependent  on  a  bulk  Richardson  number,  and  |v^|  is  the  magnitude  of 

the  mixed  layer  horizontal  velocity  (|v„j  =  '). 


2.a.(3)  Conservation  of  Potential  Temperature 


The  conservation  equation  for  moist  static  energy  is 

dh  dh  d'h  I  dQj 

dt  dx^  dxj  pc^  dx. 


(2.a.(3).l) 


where  h  =  CpT  +  gz  +  L,q  is  moist  static  energy  (L^  is  latent  heat  of  vaporization  of 
water  and  q  is  the  water  vapor  mixing  ratio),  is  the  kinematic  molecular  diffusivity  for 
moist  static  energy,  Cp  the  specific  heat  for  moist  air  at  constant  pressure,  and  Qj  the 

component  of  net  radiation  in  the  direction.  This  equation  is  then  multiplied  by  the 
mean  density  and  is  added  to  the  product  of  h  and  the  continuity  equation  which  gives 


.dp  -  dh  ,d(puj)  -  d'h 


(2.a.(3).2) 


Now,  equation  (2.a.(3).2)  is  put  into  flux  form. 


d{ph)_  d(pi4^h)  ^ 
dt  dX: 


dXj  Cp  dx^ 


(2.a.(3).3) 


Next,  equation  (2.a.(3).3)  is  expanded  into  mean  and  perturbation  parts  to  give 


^(p*) .  '^(p"/*)  p(p",*')  p(p",'*)  p(p";*') 


dx. 


dx. 


dXj  dx^ 

1  dQ^  1  dQ] 


—  d^h  _  d'h' 

dx^  dx^  Cp  dx^  Cp  dx, 


(2.a.(3).4) 


After  Reynolds  averaging,  equation  (2.a.(3).4)  becomes 
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(2.a.(3).5) 


d(ph)  ^  d(pUjh)  d(pu’h')  _  ^  I  dQj 

dt  dXj  dXj  ^  *  dx^j  Cp  dxj 


The  horizontal  derivatives  of  the  perturbation  quantities  are  neglected  using  the 
same  scaling  argument  presented  in  section  2.a.(2),  then  equation  (2.a.(3).5)  is  vertically 
integrated  ignoring  molecular  conduction  and  radiation  divergence  which  gives 


[  ^z  =  -  f 

f 

L  ^  2=1 

^  2=1 

^  2={. 

'  S{^pwh ) 


2=:s_ 


dz 


■dz 


(2.a.(3).6) 


2=2j_ 


dz 


Leibniz's  rule  is  then  applied  to  (2.a.(3).6)  to  yield 

I  lphd2  =  p*  -  p*  ^  -  -I 

OT ,  *  at  at  ax  ox 

2-!s-  2-2s- 

_ -rdZo  d  ^  Vl _ Yj _ - _ T^^S- 

puh-^-—  \pvhdz  +  pvh-^-pvh—^-  (2.a.(3).7) 

dx  dy  dy  dy 

pwh{t,  )  +  pivA  (t,  z^.)  -  pw'h'{t,  Zb^  )  + 

pw'h\t,Zs.). 


The  second  to  last  term  in  (2.a.(3).7)  is  zero  because  turbulence  goes  to  zero  above  the 
boundary  layer.  The  terms  at  the  bottom  and  top  of  the  boundary  layer  are  then  combined 
to  give 


-  lphdz  =  h 


2  =  2,. 


^  dx  dy  ^ 


Jdz 


dt 


lx 


— +  w— -^  + w 


d  e^^T. 


(2.a.(3).8) 


—  jpuhdz-—  j pvhdz  +  {^pw'h'^^  . 


2  =  2,. 


dyj 
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As  described  in  the  previous  section,  the  terms  in  the  brackets  represent  the 
entrainment  mass  flux  across  the  PEL  top  and  surface  respectively  (where  the  surface 
terms  add  to  zero  because  no  mass  can  cross  the  earth's  surface).  The  horizontal 
derivative  terms  are  combined  as  in  the  last  section,  and  then  equation  (2.a.(3).8)  reduces 
to 


3  _  _ 

—  jphdz  =  Ehg^-V»  j hpvdz  +  {^w'h'^^ 


(2.a.(3).9) 


The  integrals  are  then  evaluated  using  mixed  layer  values. 


=  Eh,,  -  V . (*,p.v.4z. )  +  [pWh'\_. 


or 


PA  ^  =  Eu„-h,  -  V  •  (*„p.v  A ) + 

(p^) 


or  employing  the  same  vector  identity  and  expression  for  the  entrainment  used  to  obtain 
equation  (2.a.(2).17), 


%  =  E(hs^  -h„)-  p„Az„v„  •  Vh„  +  (pw'h')^_. 


dt 


(2.a.(3).10) 


Now,  horizontal  homogeneity  is  assumed  in  (2.a.(3).10)  which  eliminates  the 
third  term  on  the  right  hand  side.  Then  the  entrainment  terms  are  combined,  and  the 
equation  is  divided  by  the  mixed  layer  density  and  boundary  layer  depth  which  gives 


Bh.  E(h„-h.)  y«'hX_ 


(2.a.(3).ll) 


The  equation  for  the  mixed  layer  potential  temperature  (with  no  phase  changes  of 
water)  is  based  on  the  above  equation.  It  is  obtained  by  replacing  the  moist  static  energy 
by  the  potential  temperature  and  including  the  surface  heat  flux  parameterization 
presented  in  Chapter  3  to  give 
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(2.a.(3).12) 


dQ„  ^  E{eB^-e„)  ^  p,f,a^\v„\F,^,, 

—moisture 

dt  0.74p„Az„ 


(9s--9.). 


where  is  the  same  drag  coefficient  used  in  (2.a.(2).20),  |v„|  is  defined  as  in  the 
previous  section,  and  is  a  similar  function  to 


2.a.(4)  Conservation  of  Moisture 


The  conservation  of  total  water  is  given  by 


St  ^  dXj  dxj  p 


(2.a.(4).l) 


where  is  the  molecular  diffusivity  for  water  and  5^  is  a  net  precipitation  source-sink 

term.  This  equation  is  then  multiplied  by  the  mean  density  and  is  added  to  the  product  of 
q,  and  the  continuity  equation  to  yield 


-dq,  dp  -  dq,  ^(pUj)  _ 


dt 


=  -pUj 


dXj 


dx, 


+  PK 


^-h5 

dxj 


Equation  (2.a.(4).2)  is  next  put  into  flux  form, 

-  d\  ^ 


(2.a.(4).2) 


(2.a.(4).3) 


Equation  (2.a.(4).3)  after  expanding  into  mean  and  perturbation  quantities  becomes 


d{pqf)  ^  d{pq;)  ^  djpujq,)  d{pujq;)  d(pu'jq,)  S{pUjq;)  ^ 

dt  dt  dXj  dXj  dXj  dxj 


_  d^q,  _  d^ 

1  #  *  «  I  r  - 


P^‘t. 


avi 


■p\ 


Qt 


dx^j 


+  S, 


(2.a.(4).4) 


where  the  source-sink  term,  ,  is  assumed  to  be  a  mean  forcing.  This  equation  is  then 
Reynolds  averaged  to  give 
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(2.a.(4).5) 


A  ~  axj  Ac] 

Assuming  Sg  =0  (no  precipitation  leaving  or  falling  into  an  air  parcel),  neglecting 

molecular  diffusion,  and  neglecting  horizontal  derivatives  of  perturbation  terms, 
(2.a.(4).5)  becomes,  after  vertical  integration. 


7"  -fei/z  =  -~j"  -  '7"  - 

dt  dy 


d{pwqt)^_ _ 

dz  dz 


(2.a.(4).6) 


This  equation  is  analogous  to  (2.a.(3).6).  Following  the  derivation  from  the  previous 
section,  (2.a.(4).6)  simplifies  to 


dt  P„2lz„ 


(2.a.(4).7) 


In  the  absence  of  clouds,  q,=q  (water  vapor  mixing  ratio).  Then,  equation  (2.a.(4).7)  is 
used  with  q  and  the  same  parameterization  for  the  surface  moisture  flux  as  was  used  for 
the  surface  heat  flux,  to  give 


dq„  _  ^  Pifctt  (—  _  \ 

*  ■  0.74p,az, 


(2.a.(4).8) 


2.a.(5)  Turbulence  Kinetic  Energy  (TKE)  Equation 


8- 

Equation  (2.a.(2).5)  is  expanded  into  mean  and  perturbation  parts  (except  for  the 
term  which  has  already  been  expanded  when  making  the  Boussinesq 


approximation)  to  give 


22 


(2.a.(5).l) 


dt 


dXj 


■ 

fa' V 

Si3 

g- 

ff 

. 

+  feij3(uj+u))- 


1  d{p  +  p')  ^  ^d^{Ui  +  ul) 
p  dXi  dXj 


Equation  (2.a.(5).l)  is  then  algebraically  expanded  which  yields 


V — / 

dXj 

dXj 

-^-S,3g  +  Si3 
dXj 

1  ^ 

I  ^ 

d^Ui 

d\ 

pdx, 

pdxi 

r  1 

dx] 

dx]- 

This  equation  is  then  Reynolds  averaged  and  simplifies  to 


du,  uM  u'du.  _  .  ^  1  dp  <9  M, 


av; 


(2.a.(5).2) 


(2.a.(5).3) 

Next,  the  continuity  equation  for  turbulent  fluctuations  is  multiplied  by  u'  and  Reynolds 
averaged  which  gives  m'“  =  0.  This  term  is  added  to  the  last  term  on  the  left  hand  side 

dKj 

of  (2.a.(5).3).  This  sum  is  then  put  into  flux  form, 

dil.  Ujdu-  dU-Uj)  .  _  /  ^  d% 


dt  'dx.  dx^  pdx,  dx^j 


(2.a.(5).4) 


Now,  (2.a.(5).4)  is  subtracted  from  (2.a.(5).2)  which  leaves 
du'  _  du'  ,  dUi  ,  du'  ,  c 


I  r  ,  J 


dx^j  dXj 


(2.a.(5).5) 


Equation  (2.a.(5).5)  is  next  multiplied  by  p2u'  which  gives 
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(2.a.(5).6) 


P2u'^  =  -p2uju;^  -  P2u;u'i  ^  -  p2u'l/j  ^  +  p2s„u;\ 

,,  -^u'idp'  A^j) 

P2fe,,u,u^-p2^-^p2vu.-^^p2u,—. 


1  ^ 
rJ-  g  + 

V^J 


The  Reynolds  averaged  continuity  equation,  using  mean  density,  is  then  multiplied  by 
(«,')  and  added  to  (2.a.(5).6)  which  results  in 


dXj 


dxi 


P2u’^ + wr  f -  wr  -  p2">;  ?  - 

p2u;u;  ^  +  p25uU,|&  jg  +  p2/£,jju;u}  - 


-^u[dp'  ,  A  . 

p2d:^  +  p2vu'^  +  p2u'-\^. 

p  dXi 


dxj 


dXj 


(2.a.(5).7) 


Equation  (2.a.(5).7)  is  then  put  into  flux  form. 


dt  dxj  ^  '^ar.  ^  j  ar, 

a, 


P2ulu'^-pi.j 


dXj 


p2a,in; 


a 


Ig  +  p2feijji4u)  -p2^^  +  p2vu'^  +  (2.a.(5).8) 


vy 


P  ar, 


dxj 


afu/n';) 
P2u’r  ''  " 


dXj 


Ap<).. 


The  perturbation  continuity  equation  (  —  =  0,  derived  by  using  mean  density 

av, 

in  the  continuity  equation,  expanding  this  equation  into  mean  and  turbulent  momentum 
parts,  and  subtracting  the  Reynolds  averaged  expanded  equation  from  the  expanded 
equation)  is  multiplied  by  {u'f  and  added  to  (2.a.(5).8).  The  result  is  then  put  into  flux 

form. 
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‘  ^  dXi 


St 


dXi 


p2Sisu'i 


dxj 


-.^uidp'  ,d^u' 

< _ t_ _ L  J‘\'^xrti' _ L 


p2ul 


,a{u:u'] 


dXj 


(2.a.(5).9) 


Equation  (2.a.(5).9)  is  next  Reynolds  averaged  to  give 

/p(m;)^1  ^\puj{u^^^  _ ^  ^Pu]{u^ 

^ --S:- 


Sxj 


Sxj 


p2Si3U-\  ^  |g  +  p2feij3uXj  -  +  p2vui^. 


\^vj 


(2.a.(5).10) 


Equation  (2.a.(5).10)  is  simplified  in  the  following  manner.  First,  the  second  to 
last  term  (pressure  perturbation  term)  is  rewritten  as 


(2p^d(u'p')  ^_p'Su'  ,, 

-{■—  — ^  +  2p-^— Next, 

{p  J  SXi  p  dr, 


the  repeated  indices  are  summed  over  which  eliminates  the  coriolis  term  in  (2.a.(5).10) 
and  the  last  term  in  the  expression  above  (which  converts  the  pressure  term  to  divergence 
form).  Since  TKE  is  defined  as  e  =  0.5(u'^  +  v'^  +  w'^),  it  is  appropriate  to  sum  over  the 


repeated  indices  here.  Finally,  the  last  term  (viscous  dissipation  term)  is  rewritten  as 

7 

.  The  first  term  in  this  expression  is  the  molecular  diffusion  of 

velocity  variance.  This  variance  changes  slowly  with  distance  in  the  boundary  layer  with 
typical  values  for  the  first  term  on  the  order  of  10'^^  kg  m'^  s'^.  Considering  an  eddy  0.1 
meters  in  diameter  with  a  velocity  that  changes  by  0.01  meters  per  second  across  the 
eddy,  the  instantaneous  shear  across  this  eddy  is  .1  s'^  The  shear  becomes  larger  for 
smaller  eddies.  Using  this  value,  the  second  term  is  on  the  order  of  10'^  kg  m  *  s'^.  For 
smaller  eddies  this  term  would  be  larger.  Thus,  the  first  term  in  this  expression  is  several 
orders  of  magnitude  less  than  the  second  term  and  can  be  ignored.  These  results  are  then 
applied  to  (2.a.(5).10)  which,  after  dividing  by  2,  gives 


- 


r  3..  A 


pv 


Sx^- 


-2pv 


du. 
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(2.a.(5).ll) 


d(pe)  ^(pUje)  _-r-rdUi  (uX) 

-^  = -  pu;„;  - -1-^  +  p5,,  - 


dXj  dXj 


I  dll'  I 

where  £  =  pv\  — -  .  Equation  (2.a.(5).l  1)  is  simplified  by  neglecting  horizontal 
derivatives  of  perturbation  terms.  This  gives 

d{pe)  ^  d{pue)  d{pve)  d{pwe)  d(pw'e) 


dx;  Bz 


_(w'0;)  d(wy) 


(2.a.(5).12) 


Now,  this  equation  is  integrated  from  level  S-  to  B+  (see  figure  l.b.l) 


'['■  Mla2= -'T  dz  -  T  dz  -  T  dz  - 

L  *  *  ,L  ^  .i-  * 

jpmdz-j^-ipdz-redz. 


(2.a.(5).13) 


Leibniz's  rule  is  used  to  transform  (2.a.(5).13)  to 
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-  lpedz  =  e 


dt 


I  rr  — 


dt 


■  +  u- 


dx 


+  v^^-w 


dy 


—( dz^_  ^dz^_  _dz^_  _ 

— —  +  u-z^+v — ^ — w 


dt 


dx 


dy 


— ^  I  dz-—  I  pve  dz  + 
dx .  '^  J 


Z=ZS- 


dy 

•'  2-2s- 


2  =  2b. 


-  f  pu'u'^dz+  j  p^^^gdz- 


2=2/1, 


dx. 


2=rj_ 

2=2*, 


e.. 


|-  djpw'e)  -i-  a{wy) 


2=2*, 


2  =  2*- 


2=2*- 


dz 


dz-  J e dz. 


2  =  2*- 


(2.a.(5).14) 


This  equation  is  similar  to  (2.a.(2).14),  where  the  bracketed  portion  of  the  first 
two  terms  on  the  right  hand  side  is  the  entrainment  mass  flux  across  the  B+  and  S- 
surfaces.  The  first  term  is  zero  because  the  turbulence  kinetic  energy  vanishes  just  above 
the  PBL  top,  and  the  second  term  is  zero  because  the  entrainment  mass  flux  across  the 
earth's  surface  is  zero.  The  seventh  term  on  the  right  hand  side,  which  represents  the  flux 
divergence  of  TKE,  is  zero  because  the  vertical  turbulent  flux  of  TKE  is  equal  at  the  top 
and  bottom  of  the  PBL,  hence  the  vertical  integration  of  this  quantity  is  zero.  Finally,  the 
second  to  last  term  on  the  right  hand  side  (pressure  correlation  term)  is  zero  since 
turbulence  vanishes  at  the  surface,  and  if  gravity  waves  (that  remove  TKE  from  the  top  of 
the  PBL)  are  neglected.  Equation  (2.a.(5).14)  is  then  rewritten,  using  the  hydrostatic 

relation  (^  =  -pg), 
dz 


2  =  2„ 


^  jpedz  =  -V*  jepvdz-  J  F^*^dz-  J  F^»^£/z  + 


dv 


dt 


2  =  2j_ 


2  =  2*- 


dx 


dy 


P=Ps- 


I  J 


p-  — 

P=Pb.  V 


P=Pb* 


dp 


J  2-2*.^ 

-dp-  \edz, 

P  2  =  2, 


(2.a.(5).15) 


where  =  pv'u',  F^  =  pv'v',  and  F^  =  pv'w'  are  each  three  components  of  the  turbulent 
momentum  flux  (also  known  as  the  Reynolds  stress).  The  subscript  v  (for  vector)  is  used 
instead  of  w  in  the  last  equality  because  this  quantity  will  be  called  the  vector  momentum 
flux.  The  first  and  second  integrals  are  simplified  by  using  the  mixed  layer  values  for  the 
density,  TKE,  and  momentum  which  gives 
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or 


^  m/  J  *  x 

Z=2s~ 

dv 


2=2b* 


J 


2s-  P=P». 


P=Pb* 


^  d. 


>\\  J  2=^0. 

—dp-  [edz, 
P  ^ 

^  2=Zj- 


_  [  f  f  f  dz 

J^;iy  J  >  A, 


2  =  2s, 

P=Ps 


dx 


dy 


P=Ps-  f  P=Ps- 

1  I 

P=Pb+  ^  ^  P=Pb* 


2=2S- 

P=Ps-/'  f,..'a>\\ 


J'-y:) 

^  0.. 


2  =  28* 


-dp-  \edz. 

P  ■' 


Z=2s- 


(2.a.(5).16) 


The  shear  term  is  then  simplified  by  breaking  the  integration  up  into  surface 
layer,  mixed  layer,  and  entrainment  layer  components.  For  the  surface  layer,  the  flux  is 
constant  and  parallel  to  the  wind  (either  both  quantities  are  negative  or  positive,  therefore, 
the  product  and  integration  of  the  product  are  positive).  The  wind  increases  from  zero  at 
S-  to  its  mixed  layer  value  at  S.  Therefore, 


P-Ps  11— I 


^Top  of  Surface  Layer 


P-Ps- 


or 

^Top  of  Surface  Layer  |(^ 


(2.a.(5).17) 


The  negative  sign  in  the  first  equation  of  (2.a.(5).17)  is  needed  because  the  limits  of 
integration  were  reversed,  but  as  was  stated  above,  the  result  (second  equation)  is 
positive.  For  the  mixed  layer,  the  wind  is  equal  to  its  mixed  layer  value,  thus  there  is  no 
shear  here.  In  the  entrainment  layer,  the  flux  is  assumed  to  decrease  linearly  from  its 
value  at  B  to  zero  at  B+.  The  wind  changes  from  its  mixed  layer  value  at  B  to  another 
value  at  B+.  Again,  the  wind  and  flux  are  assumed  to  be  parallel.  Then, 
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^Top  of  Entrainment  Layer  2 

or 


P=Pb 


^Top  of  Entrainment  Layer 


(2.a.(5).18) 


where  |4v|  =  |va+|-lv„|  and  (f,)^  =  -E\Av\.  The  second  equality  in  (2.a.(5).18)  is  valid 

because  the  flux  of  momentum  into  the  top  of  the  boundary  layer  is  due  to  entrainment  of 
air  from  the  free  atmosphere  when  there  is  wind  shear  through  the  entrainment  layer. 

This  flux  is  zero  if  the  entrainment  rate  is  less  than  zero  or  if  there  is  no  wind  shear.  The 
second  equality  is  obtained  in  the  following  manner.  Neglecting  fluxes  due  to  radiation 
and  clouds,  the  rate  at  which  mass  is  added  to  the  PBL  from  the  free  atmosphere  (FA)  is 
given  by  gE.  For  any  arbritary  variable.  A,  the  upward  turbulent  flux  of  A  is  denoted  by 
Fa-  The  continuity  of  the  total  flux  at  level  B,  assuming  =  0,  is 

-EAg_^  =  -EAg  +  {Ffjg.  The  flux  added  to  the  PBL  from  the  free  atmosphere  must  equal 

the  total  flux  within  the  PBL  which  consists  of  the  flux  within  the  PBL  due  to  mass 
entrainment  and  the  flux  within  the  PBL  due  to  upward  turbulence  transport.  The 
transition  of  A  at  the  PBL  top  is  modeled  as  a  jump  given  by  A4  =  Ag^  -  Ag.  Using  this 
in  the  flux  continuity  equation  gives  {Ff)g  =  -EM.  If  A  =  v  then  the  second  equality  is 

obtained.  Next,  horizontal  homogeneity  is  assumed  in  equation  (2.a.(5).16),  except  for 
the  second  term  which  contains  the  mean  divergence.  This  eliminates  the  third,  fourth, 
and  fifth  terms  in  this  equation.  Then,  (2.a.(5).17)  and  (2.a.(5).18)  are  summed  with  the 
result  substituted  into  (2.a.(5).16)  to  give 


=  -«.P.4z.(  V  .  V.) + 1  £|  Af  +  |(F.  ),|v., 


P-Ps- f 


I 


P-Pb* 


pm 

e.. 


Z  =  Z«. 


—dp-  {edz. 
P  ^ 


Z=Zs- 


(2.a.(5).19) 


The  buoyancy  term  can  also  be  simplified  by  using  an  approximation  for  the  flux 
of  virtual  dry  static  energy.  The  buoyancy  term  is  written 


P=Ps- 


B=  I 


P=Pb. 


om'' 

^  e.. 


1  . 
—dp, 

P 
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or  using  Fg  =  pw'd' , 


P=Ps- 


B  =  J^ldp  =  -L  J 


P^Ps* 


P  P^Pb* 


fpZA 

B^ 


where 

w 


or  using  = 


and  equation  of  state. 


B  = 


P-Pb- 


-  I 

c„  •' 

P  P=Pb, 


J7  ^  17 

^dp  =  K  f  ^4?. 

P  •’  P 


(2.a.(5).20) 


Since  the  pressure  does  not  change  that  much  in  the  PBL  and  the  flux  of  virtual  dry  static 
energy  varies  linearly,  their  ratio  is  nearly  linear,  and  the  last  integral  in  C2.a.(5).20)  can 
be  simplified  to 


i^Sv)s  j  (^5v)fl 

Ps  Pb 


(2.a.(5).21) 


Now,  the  definition  of  the  flux  of  virtual  dry  static  energy  from  (2.a.(5).20)  is  used  in 
(2.a.(5).21)  to  get 


'cAP.)sf 

/  \ 
Pe 

\ 

PsB„ 

H 

or  using  Poisson's  equation. 


5  =  ^ 
2 


_L. 

(Pb) 

[[UJ  PsB^ 

[Po) 

PB^m  J 

or 


2 


'ftTM 

^Po)  Ps 


Po 


Pb 
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or  using  =  F^  +  {e„\5F^,  where  ^  =  -~  =  ^  'q^j2  ^  ’ 


UJ  Pb 


(2.a.(5).22) 


where  (Fg)^,  (Fg)^  =  -EAd,  (F^)^,  and  (F^)^  =  -EAq  are  the  PBL  surface  and  PBL  top 

heat  and  moisture  fluxes  respectively.  The  PBL  top  fluxes  of  heat  and  moisture  are 
defined  in  the  same  manner  as  the  mass  flux  into  the  PBL  top.  Chapter  3  contains  the 
details  of  the  surface  heat  and  moisture  flux  parameterizations. 

The  dissipation  term  is  modeled  by  using  a  second-order  closure  assumption. 

2=2e, 

The  vertically  integrated  dissipation  rate  is  D  =  jedz  =  p„(j\  where  a  is  the  dissipation 

velocity.  Closure  is  obtained  by  assuming  the  square  of  o  is  proportional  to  the  vertically 
averaged  TKE,  <j'  =  —,  where  a,  =  0.163  based  on  Deardorff’s  (1974)  results.  The 

a, 

vertically  integrated  dissipation  rate  is  then  related  to  the  vertically  averaged  TKE  by 


D  =  p„  — 

a, 


(2.a.(5).23) 


The  TKE  equation  is  finally  wiiaen 


d{p„e„Az^) 


=  •  v„)  +  S-i-  5  -  D, 


Az  p  ^^  =  -e  r +  p  ^  (V*v  )  +S  +  B-D, 

mr  m  ^  ^  mf  ^ 
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or  using  the  hydrostatic  relation  and  since  the  term  inside  the  brackets  is  the  entrainment. 


%  =  -^(S  +  5  -  -  D),  (2.a.(5).24) 

where  S  is  given  by  (2.a.(5).17)  and  (2.a.(5).18),  B  by  (2.a.(5).22),  and  D  by  (2.a.(5).23). 
The  entrainment  rate  is  determined  in  the  following  section. 

2.a.(6)  Entrainment  Rate  Equations 

Entrainment  is  the  mechanism  that  brings  unmixed  free-atmosphere  air  into  the 
top  of  the  PBL.  This  air  becomes  mixed  by  the  existing  turbulence  in  the  mixed  layer 
causing  the  mixed  layer  to  grow.  The  entrainment  rate  is  positive  if  free-atmosphere  air 
is  being  brought  into  th'"  top  of  the  mixed  layer.  It  is  zero  if  no  air  is  transported  across 
the  PBL  top.  If  air  is  being  removed  from  the  top  of  the  mixed  layer  then  the  entrainment 
rate  is  negative  and  the  mixed  layer  is  decaying.  Because  the  previously  described 
prognostic  equations  require  knowledge  of  E,  it  must  be  parameterized  to  solve  these 
equations.  Since  turbulence  is  required  to  mix  newly  entrained  free-atmosphere  air,  E  is 
considered  proportional  to  the  square  root  of  the  TKE.  This  is  the  basis  for  the 
parameterization  described  in  section  2.a.(6)(a).  This  parameterization  is  used  when  the 
entrainment  rate  is  determined  to  be  positive. 

If  there  is  no  turbulence,  then  the  TKE  and  E  will  equal  zero.  The  existence  of 
turbulence  alone,  however,  does  not  guarentee  that  E  will  be  positive.  Table  2.a.(6).l 
summarizes  the  conditions  that  determine  the  sign  of  E. 
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Table  2.a.(6).l;  Sign  of  E  Based  on  B+S  and  Bq+Sq- 


Sum  of  Buoyancy  and  Shear  Computed  with  E=0 
(where  0  <  fraction  <  1) 

Sum  of  Buoyancy  and 

Shear 

Bfl  +  So  <  fraction  *  D 

Bo  +  So  >  fraction  *  D 

Case  1 

Case  2 

E<0 

E>0 

B  +  S  <  fraction  *  D 

(B  <  Bo  /.  EB,  <  0;  since 
EB,  =  0  if  E  <  0,  E  must 
be  greater  than  zero  and 

B,  must  be  less  than  zero) 

Case  3 

Case  4 

E>0 

E>0 

B  +  S  >  fraction  *  D 

(B  >  Bo  A  EB,  >  0;  since 
EB,  =  0  if  E  <  0,  E  must 
be  greater  than  zero  and 
Bj  must  be  greater  than 
zero) 

These  conditions  are  checked  during  each  time  step  of  a  model  run.  If  E  is  determined  to 
be  less  than  zero  during  a  given  time  step,  then  the  negative  parameterization  described  in 
section  2.a.(6)(b)  is  used  to  compute  E. 

During  rapid  growth,  the  TKE  and  E  become  large.  Since  the  dissipation  rate  is 
proportional  to  the  TKE,  it  also  becomes  large.  It  is  possible  for  the  sums  of  B+S  and 
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Bq+So  to  be  less  than  D  which  would  cause  a  large  entrainment  rate  to  suddenly  become 
negative.  There  are  two  ways  to  prevent  this  from  happening.  One  way  is  to  set  a 
threshold  of  E  such  that  when  E  is  currently  greater  than  this  threshold  it  is  calculated 
using  the  positive  parameterization  regardless  of  the  value  of  the  sums.  A  better  way  is  to 
check  the  sums  against  a  fraction  of  the  dissipation  (as  indicated  in  Table  2.a.(6).  1), 
where  this  fracuon  is  set  as  a  tunable  parameter.  A  proper  selection  of  this  parameter  will 
prevent  the  sums  from  being  less  than  the  fraction  of  the  dissipation  during  periods  of 
rapid  PBL  growth.  Table  2.a.(6).l  is  discussed  further  in  section  2.a.(6)(b). 

During  the  late  afternoon,  before  sunset,  the  clear  convective  boundary  layer  over 
land  has  reached  a  quasi-steady-state.  At  this  point  the  surface  buoyancy  flux  rapidly 
approaches  zero  with  the  loss  of  daytime  heating.  Both  the  entrainment  rate  and  TKE  are 
small  compared  to  their  values  in  the  mid-morning  (during  the  rapid  growth  of  the  PBL). 
At  this  point,  a  balance  has  occurred  in  the  TKE  equation  Since  there  are  no  processes  to 
generate  a  significant  amount  of  TKE  at  this  time  of  day,  the  local  rate  of  change  of  the 
TKE  is  small  and  can  be  neglected.  The  sign  of  the  entrainment  rate  then  depends  on  the 
sum  of  the  buoyancy  and  shear  terms.  If  this  sum  is  small  enough  the  entrainment  rate 
will  be  negative. 

There  are  contributions  to  the  buoyancy  and  shear  production  from  the  surface 
and  PBL  top.  The  contributions  from  the  PBL  top  depend  on  the  sign  and  magnitude  of 
the  entrainment  rate.  Since  the  entrainment  rate  determines  how  fast  mass  is  brought  into 
the  top  of  the  PBL,  mass  will  cross  the  free-atmosphere  PBL  top  interface  only  when  the 
entrainment  rate  is  positive.  Thus,  if  the  entrainment  is  zero  or  less  then  there  will  be  no 
contribution  to  the  buoyancy  or  shear  production  at  the  top  of  the  PBL.  With  positive 
entrainment,  buoyancy  production  at  the  top  of  the  PBL  can  be  positive  or  negative 
depending  on  the  gradient  of  temperature  and  moisture  here.  The  shear  production  at  the 
PBL  top  is  always  non-negative.  It  is  positive  if  the  entrainment  rate  is  positive  and  there 
is  wind  shear  across  the  top  of  the  PBL,  and  it  is  zero  if  either  the  wind  shear  is  zero  or 
the  entrainment  rate  is  zero  or  less.  Therefore,  the  sum  of  the  buoyancy  and  shear 
production  (with  surface  and  top  contributions)  along  with  the  sum  of  the  buoyancy  and 
shear  production  at  the  top  of  the  PBL  must  be  considered  to  determine  the  sign  of  the 
entrainment  rate. 
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2.a.(6)(a)  Positive  Entrainment 


From  Breidenthal  and  Baker  (1985),  the  positive  entrainment  rate  formula 
without  clouds  is 


bi 

{l  +  b2Riy 


(2.a.(6)(a).l) 


where  Ri  =  js  relevant  Richardson  number,  and  bi 

determined  as  follows.  For  a  strong  inversion,  ^2/?/ » 1,  and  (2 


Now,  substituting  the  expression  for  Ri  into  (2.a.(6)(a).2)  gives 
b,  _  gEAd^Az^ 


and  b2  are  constants 
a.(6)(a).l)  reduces  to 

(2.a.(6)(a).2) 


(2.a.(6)(a).3) 


Using  the  famous  “0.2”  formula,  EA9^  =  0.2(Fe^  (see  Randall,  1984),  (2.a.(6)(a).3) 
becomes 


b,  _  0-2g(F9j/z„ 


(2.a.(6)(a).4) 


Next,  a  balance  is  assumed  between  buoyant  production  and  dissipation  of  TKE.  The 
buoyancy  term  is  written  in  a  slightly  different  form  from  (2.a.(5).22),  and  the  dissipation 
is  given  by  (2.a.(5).23).  This  balance  is  then  written 


(PeJ 

1  Az^  ; 

c  1 

203  ‘■^"’1 

Ui  J 

(2.a.(6)(a).5) 


Since  =  Pn,  and  (0„)^  =  0s,  (2.a.(6)(a).5)  can  be  substituted  into  (2.a.(6)(a).4)  to 
obtain 
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(2.a.(6)(a).6) 


bi  _  {02)2 


6.10. 


Now,  in  the  no  inversion  limit  (Ri=0),  (2.a.(6)(a).l)  reduces  to 
E  =  PB4^xbi- 


(2.a.(6)(a).7) 


Deardorff  (1974)  found  by  large-eddy  simulation  that 


E  =  0.2pB 


(p8)s 


;1/J 


(2.a.(6)(a).8) 


where  w 


(P9), 


ni/3 


is  the  convective  velocity  scale  of  Deardorff  (1970).  Using 


(2.a.(6)(a).5)  and  (2.a.(6)(a).8)  (where  Ps=  Pm=  Pb)>  (2.a.(6)(a).7)  becomes 


bi  =  0.2 


(«/) 


3/2 


«  0.624. 


(2.a.(6)(a).9) 


Finally,  from  (2.a.(6)(a).6),  one  obtains  62  ~  0.102 . 


2.a.(6)(b)  Negative  Entrainment 

Assuming  the  entrainment  rate  and  TKE  are  small  compared  to  their  values 
during  rapid  PEL  growth,  the  local  rate  of  change  of  TKE  is  small  and  can  be  neglected 
in  the  TKE  equation  to  give 

0  =  -Ee„  +  S  +  B-D.  (2.a.(6)(b).l) 


The  sign  of  the  entrainment  rate  then  depends  on  the  sum  of  the  buoyancy  and  shear 
terms.  Solving  for  E  in  the  above  equation  gives 


£  = 


S  +  B-D 


(2.a.(6)(b).2) 
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where  E>0  if  (S+B)>D.  There  are  four  possible  cases  for  determining  the  sign  of  E.  The 
buoyancy  and  shear  terms  are  first  written 

B  =  Bo  +  EB,  and  S  =  So  +  £S, ,  (2.a.(6)(b).3) 

where  the  zero  subscript  indicates  the  surface  contribution  to  the  buoyancy  and  shear  (as 
if  E  were  zero),  and  the  one  subscript  is  the  contribution  to  the  buoyancy  and  shear  at  the 
top  of  the  PBL  due  to  entrainment  (as  if  the  surface  fluxes  were  zero).  The  buoyancy 
and  shear  terms  (B  and  S)  are  then  computed  assuming  E>0.  Then,  the  surface 
contributions  to  the  buoyancy  and  shear  (Bq  and  So)  are  computed  and  summed.  These 
sums  are  compared  to  arrive  at  one  of  the  four  possible  cases  listed  in  Table  2.a.(6).l. 

If  Case  1  occurs  then  the  entrainment  rate  is  determined  using  the  negative 
production  formulation.  This  is  accomplished  by  partitioning  the  TKE  equation  into  a 
weighted  contribution  of  the  local  rate  of  change  of  TKE  and  a  weighted  contribution  of 
the  production  of  TKE  due  to  entrainment.  Equation  (2.a.(5).24)  is  split  into  two 
equations  (where  B=Bo  and  S=So  since  E<0), 

^  vvcig/it(Bo  +So-D) 

and  (2.a.(6)(b).4) 

^  (l-wcig/it)(Bo+5o-D) 

t,  —  - , 

where  0  <  weight  ^  1 .  If  the  weight  is  set  to  one  then  the  sum  of  the  above  equations  is 
just  (2.a.(5).24).  The  TKE  is  first  determined  using  the  top  equation  in  (2.a.(6)(b).4),  and 
then  the  entrainment  rate  is  determined  using  this  new  value  of  the  TKE  and  the  bottom 
equation  in  (2.a.(6)(b).4). 

2.b.  Initialization 

The  model  requires  that  certain  variables,  including  prognostic  variables,  be 
initialized  before  time-stepped  predictions  are  made.  Chapter  4  provides  a  brief 
description  of  the  Wangara  data  set  used  to  initialize  the  land  simulations.  For  land,  four 
data  files  are  used  that  include:  three  hourly  sounding  data  which  includes  temperatures 
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and  mixing  ratios  at  various  pressure  levels,  hourly  sounding  data  which  includes  u  and  v 
wind  components  at  various  heights,  hourly  ground  temperatures,  and  hourly  u  and  v 
geostrophic  wind  components.  Prognostic  variables  initialized  over  land  are  obtained  by 
interpolating  between  two  data  periods  based  on  the  model  start  time  (e.g.,  with  a  start  of 
1030L,  temperatures  and  mixing  ratios  would  equal  the  sum  of  one-half  of  their  values  at 
the  0900L  and  1200L  sounding  times),  and  by  interpolating  between  data  levels  (heights 
or  pressures)  where  appropriate.  Table  2.b.l  sununarizes  the  prognostic  variables  that  are 
initialized  for  simulations  over  land  or  water. 


Table  2.b.l:  Summary  of  Prognostic  Variable  Initializations. 


Prognostic  Variable 

Land  Initialization 

Water  Initialization 

Mean  pressure 

thickness 

Initialized  based  on 

data  and  start  time 

of  model  run 

Assigned  an  initial 

value 

Mean  u  and  v  wind 

components 

Initialized  from  data 

Computed  as  one- 

half  sum  of  surface 

wind  and  wind  at 

top  of  PEL 

Mean  potential 

temperature 

Initialized  from  data 

Assigned  an  initial 

value 

Mean  mixing  ratio 

Initialized  from  data 

Assigned  an  initial 

value 

Turbulence  kinetic 

energy 

Assigned  an  initial 

value 

Assigned  an  initial 

value 

Table  2.b.2  contains  the  constants  that  must  be  set  at  the  start  of  a  simulation. 
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Table  2.b.2:  Summary  of  Constants. 


Constant 

Land  Value 

Water  Value 

Use 

Time  step  (At) 

60  seconds 

60  seconds 

Compute  prognostic 
equations 

Sea  surface 

temperature 

N/A 

289°  K 

Compute  surface 
potential 

temperature,  surface 
mixing  ratio,  and 
surface  density 

Mixing  ratio  at  top 

of  PBL 

Not  a  constant: 

interpolated  between 
two  sounding 
periods 

1  g/kg 

Compute  virtual 
potential 

temperature  at  top  of 
PBL,  buoyancy,  and 
mean  mixing  ratio 

Potential 

temperature  lapse 

rate  above  the  PBL 

N/A 

4°  K/km 

Compute  potential 
temperature  at  top  of 

PBL 

Surface  u  and  v 

wind  components 

N/A 

2  m/s 

Compute  u  and  v 
wind  components  at 
top  of  PBL 

Wind  lapse  rate 

above  the  PBL 

N/A 

5  m/s/km 

Compute  u  and  v 
wind  components  at 
top  of  PBL 

Geostrophic  u  and  v 
wind  components 

Not  constant: 

interpolated  between 
two  sounding 
periods 

Ug=-10  m/s 

Vg=0  m/s 

Compute  mean  u 

and  V  wind 

components 

2.C.  Top  Boundary  Conditions 

As  shown  by  figure  l.b.l,  the  model  domain  is  bounded  at  the  top  by  the  free 
atmosphere  and  at  the  bottom  by  the  earth's  surface.  Top  boundary  conditions  are  applied 
at  the  PBL  top-free  atmosphere  interface,  and  surface  boundary  conditions  are  applied  at 


the  PBL  bottom-earth  surface  interface.  Lateral  boundary  conditions  are  not  required 
with  the  assumption  of  horizontal  homogeneity. 

Turbulence  in  the  mixed  layer  results  in  uniform  prognostic  variables  within  the 
layer  (i.e.,  d  ^9^,  e—^e^,  «  ->  u^,  v  ->  ).  This  turbulence  also  mixes  free- 

atmosphere  air  that  is  entrained  into  the  top  of  the  PBL.  The  first  and  very  important  top 
boundary  condition  that  the  model  requires  is  that  the  turbulence  becomes  zero  at  the 
interface  between  the  PBL  top  and  the  free  atmosphere.  This  boundary  condition  is  used 
to  simplify  the  prognostic  equations.  Zero  flux  at  the  interface  is  the  boundary  condition 
that  leads  to  =  -£AA  (described  in  section  2.a.(5)).  The  next  boundary  condition 

is  applied  to  the  flux  of  A  at  level  B,  not  at  the  interface,  when  the  entrainment  rate  is  less 
than  zero.  This  flux  is  zero  when  the  entrainment  rate  is  less  than  zero  because  no  mass 
enters  the  PBL  top  when  E<0.  The  model  uses  this  boundary  condition  to  set  the  fluxes 
of  heat,  moisture,  and  momentum  across  the  PBL  top  to  zero  whenever  the  entrainment  is 
less  than  zero.  The  remaining  top  boundary  conditions  are  applied  to  the  TKE  equation. 
Since  TKE  is  a  measure  of  the  turbulence,  the  TKE  also  vanishes  at  the  interface.  The 
next  boundary  condition  applies  to  both  the  top  and  bottom.  The  vertical  turbulent  flux 
of  TKE  at  the  top  is  equal  to  its  value  at  the  surface.  The  final  top  boundary  condition  is 
that  the  pressure  correlation  term  vanishes  at  the  top  of  the  PBL  when  gravity  waves  are 
neglected. 

2.d.  Surface  Boundary  Conditions 

The  earth's  surface  acts  as  physical  barrier  at  the  bottom  of  the  PBL.  Complications 
arise  in  applying  surface  boundary  conditions  when  the  surface  is  heterogeneous  and 
varies  orographically.  The  first  surface  boundary  condition  is  horizontal  homogeneity. 
The  next  boundary  condition  is  that  no  mass  can  cross  the  earth's  surface.  These  two 
boundary  conditions  are  used  with  all  the  prognostic  equations.  The  equality  of  the 
vertical  turbulent  flux  of  TKE  at  the  top  and  bottom  is  the  third  surface  boundary 
condition.  The  remaining  boundary  condition  is  the  loss  of  turbulence  at  the  earth's 
surface.  This,  along  with  the  neglect  of  gravity  waves  mentioned  above,  eliminates  the 
pressure  correlation  term  in  the  TKE  equation.  Table  2.a.l  lists  the  boundary  conditions 
used  with  the  prognostic  equations. 
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3.  One-Layer  Model  Time  Schemes 


3.a  Surface  Heat-Moisture  and  Momentum  Flux  Parameterizations 

The  parameterization  scheme  developed  by  Louis  (1979)  is  well  suited  to  this  model. 
The  parameterization  is  complicated  enough  to  accurately  represent  the  effects  of 
boundary  layer  fluxes  over  long  periods,  but  not  too  complicated  to  preclude  rapid 
computer  solutions  even  with  lengthy  simulations.  This  is  especially  important  in 
incorporating  this  model  into  a  general  circulation  model  where  very  long  simulation 
periods  are  required.  The  parameterization  also  fits  well  with  the  boundary  layer  being 
represented  by  one  or  two  levels,  and  the  assumption  that  the  fluxes  vary  linearly  with 
height  from  the  surface  to  the  PEL  top  (the  top  may  be  constant,  or  in  this  case 
prognostically  determined  by  the  model).  The  description  of  the  boundary  layer  by  this 
model  is  sufficiently  detailed  to  prevent  incorrect  feed-back  from  occurring.  Accurate 
feed-back  is  necessary  because  this  parameterization  relates  the  magnitude  of  the  fluxes 
to  the  prognostic  variables.  Finally,  since  the  model  depends  on  both  buoyant  and  shear 
driven  turbulence,  the  parameterization  should  simulate  both  of  these  processes.  Louis' 
parameterization  accomplishes  this  by  requiring  that  the  diffusion  coefficients  not  only 
depend  on  the  wind  shear,  but  also  the  static  stability  of  the  atmosphere. 

The  parameterization  scheme  is  based  on  Monin-Obukhov  similarity  theory.  The 
Monin-Obukhov  scale  height  is  given  by 


where  is  the  scaling  velocity,  k  the  Von  Karman  constant,  g  the  acceleration 

of  gravity,  and  0,  =  -w'd'/u.  the  scaling  temperature.  The  integrated  flux  profile 
relationships  give 

u  =  j[ln(z/z„}-  (z/L)  +  ( z^  /  L)] ,  (3. a. 2) 
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and 


Q 

=  -R-j^{ln(zlZ„)  -  W,eat-nu.is,uMIL}  +  ¥ ,ea,-^is,ure( ^0 


(3.a.3) 


where  z  is  the  surface  height  (set  to  10  meters  in  the  model),  zq  the  roughness  height,  \|/ 
are  Businger's  functions  for  momentum  and  heat-moisture,  R  is  a  constant  equal  to  0.74, 
and  Ad  =  d^_  -  0^  (the  opposite  of  Louis'  definition,  hence  the  minus  sign  in  (3.a.3)). 

Substitution  of  (3.a.2)  and  (3.a.3)  into  (3.a.l)  results  in 


du^  [ln(zlzJ-\lf^JzlL)+\ifHea,-n,ois,ure(^0  '  L)] 
gAd  [ln(  zizo)-  (z/L)  +  ( Zg  /  L)f 


(3.a.4) 


The  momentum  and  heat  flux  formulations  are  then  determined  in  the  following 
manner.  First,  the  square  of  the  scaling  velocity  is  solved  for  from  (3. a.  1).  This  gives 


2  _  kgd,L 

a 


(3.a.5) 


Next,  (3.a.4)  is  substituting  into  (3.a.5)  to  get 


^  ^  ku^d.  [ln(z/Zg)-ii/„„Jz/L)+y/,^^,_„^i^,^JZo/L)] 
Ad  [ln(zlZo)-\{/^JzlL)  +  \f/„^JZglL)]' 


(3.a.6) 


Then,  equation  (3.a.3)  is  used  to  convert  (3.a.6)  to 


(3.a.7) 


where  a  = 


[in(zlzo)] 


is  the  drag  coefficient,  and  F  a  function  dependent  on  z,  zq,  and  L. 


The  Monin-Obukhov  scale  height  is  related  to  the  bulk  Richardson  number. 

Rig  =  -  ^ ,  which  can  be  inferred  from  (3.a.4),  therefore  F  is  also  dependent  on  z,  zq, 

du 

and  RL.  Thus,  the  surface  momentum  flux  is 


Ps-w'u'  =  -ps.id  =  -ps-a'u^F^„{zlZg,Rig). 


(3.a.8) 


Similarly,  the  surface  heat  and  moisture  fluxes  are 
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(3.a.9a) 


Ps-w'6'  =  -psju^dt  = 


Ps-a'uAdF^^,_^.^,^^{z  /  Zq  ,  /?/fl ) 


R 


and 


-7-,_  ^  _  Ps-a'uAqFHea,-^is,ure  '  ^0’^^  ) 


Ps-^Q  =-ps-^*Q*  = 


R 


(3.a.9b) 


The  model  uses  ( )  =  I v„  |m„  ,  ) 

V  ! u-momentum  I  m|  m’  V  / 

the  above  equations. 


V— momentum 


=  l»’mk,and  (m  ) 


= 

heat— moisture  ' 


Louis  computed  the  momentum  and  heat-moisture  functions  numerically. 
Analytical  formulae  were  then  fit  to  the  functions.  The  analytical  formulae  eliminate  the 
need  to  perform  an  iterative  calculation  during  each  time  step.  For  unstable  conditions 
(when  Rig  <0) 


F  =  l- 


bRL 


1  +  clRig 


\II2  ’ 


(3.a.l0) 


where  b=9.4,  =  7.4a^b\  —  ,  and  =  5.3a^b 


y/2 
o) 


—  .  The  function  for  the 

\^o) 


neutral  and  stable  cases  {Rig  >  O')  is 


f’  = 


U  +  b'Rigf’ 


(3.a.ll) 


where  b'  =  4.7. 

3.b.  Conservation  of  Momentum 


The  conservation  of  momentum  equations  (2.a.(2).20)  and  (2.a.(2).21)  are 
approximated  by  a  forward  time  scheme  for  the  first,  101st,  201st,  301st,  etc.,  time  steps, 
and  by  a  leap-fi’og  time  scheme  for  all  other  time  steps.  This  is  illustrated  in  figure  3.b.  1 
below.  Periodically  using  a  forward  time  step  prevents  any  large  divergence  from 
building  up  in  the  solutions  produced  by  the  leap-frog  time  steps. 
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F:  Forward  Time  Step 
L:  Leap-Frog  Time  Step 


FI  L2  L4...  LlOO  FlOl  L102..L... 


TTTYm 


starts  over  here 


Figure  3.b.  1 :  Conservation  of  Momentum  Time  Scheme  Flow  Diagram 


The  forward  time  difference  schemes  for  (2.a.(2).20)  and  (2.a.(2).21)  when  E>0  are 


I  -  W  p  ^  p  , 

rm^^m  rm^^m  ' 


or 


c'+4(f(vr-v,). 


AtEu, 


B+ 


U„  = 


aA 


"  ^  C‘  ,  AtE  ' 

I  nt  nt  t  rn  tn 


(3.b.l) 


and 


Aip 


'd . 


■  + 


p  Az 


or 
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^  ^  ^  AtE 


(3.b.2) 


For  E<0,  equations  0.b.l)  and  (3.b.2)  reduce  to 


c'+Mc'-v,) 


(3.b.la; 


P  42 

r  m  m 


r3.b.2a) 


The  pressure  gradient/coriolis  terms  are  represented  explicitly  (n-1),  while  the  divergence 
and  flux  terms  are  represented  implicitly  (n).  Fully  implicit  representation  would  require 
solving  two  equations  in  two  unknowns  simultaneously.  The  partially  implicit 
representation  is  used  to  simplify  solving  the  equations  and  still  maintain  stability.  Initial 
condition  data  (see  section  2.b.)  is  used  for  the  values  at  n-1  for  the  first  time  step.  The 
values  computed  by  the  previous  leap-frog  time  step  are  used  as  initial  conditions  (n-1 
values)  for  the  forward  time  step  computations  at  time  steps  101, 201,  301,  etc. 

The  leap-frog  time  difference  schemes  for  (2.a.(2).20)  and  (2.a.(2).21)  when  E>0  are 
=  uZ‘  +  2Atf(v:  -v\-  I*'"’  I  +  -C) 
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For  E<0,  0.b.3)  and  (3.b.4)  become 


and 


«r'  = 


iC'+w(v:-v,) 


7  + 


2AtPsf,a^F^„  V, 


vr'=- 


v;-'-2/iff;  (»;-»,)+ 


2Ata'F^^  v„ 

mom  m 


-/ 


P.^2a 


^  ^  22itp^c^'F„, 


.n-J 


Pm^n 


(3.b.3a) 


(3.b.4a) 


Here  again,  the  pressure  gradient/coriolis  terms  are  represented  explicitly  (n),  and  the 
divergence  and  flux  terms  are  represented  implicitly  (n+1).  After  the  values  at  n+1  are 
computed  during  a  leap-frog  time  step,  the  values  at  n-1  are  updated  to  the  values  at  n  and 
the  values  at  n  are  updated  to  the  newly  computed  values  at  n-i-1,  before  the  next  leap-frog 
time  step. 


3.C.  Conservation  of  Potential  Temperature 


Equation  (2.a.(3).12) 


~  moisture 

^  0.74P„AZ„ 


rewritten  as 


(3.C.1) 
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where  V  =  jg  ventilation  or  surface  layer  mass  flux.  A  backward 

0.74 

(implicit)  time  scheme  is  used  to  represent  (3.C.1).  This  scheme  is  unconditionally  stable 
and  has  first  order  accuracy.  The  finite  difference  form  of  (3.c.  1)  when  E>0  is  then 


^n+/  __  0n  gEAlf^ 


B+ 


-0^(9.- -er'). 


or 


6"^'  = 


Ap„ 


Ap„, 


s- 


J  I  gEAt  ^  gVAt 


Ap„  Ap„ 


(3.C.2) 


When  E<0  0.C.2)  reduces  to 

Qn  gV^Q 

Ap„ 

Apn, 


(3.c.2a) 


Equations  (3.C.2)  and  (3.c.2a)  were  used  for  the  ocean  simulations. 

A  prescribed  surface  heat  flux  was  used  for  the  Wangara  simulations.  Following  Andre 
et  al.  (1978),  the  surface  heat  flux  is  approximated  from  Day  33  of  the  Wangara  data  as  a 
sine  wave  (Figure  3.c.  1) 


600 


(3.C.3) 


where  is  the  kinematic  heat  flux  in  units  of  K  m  s  *,  the  maximum  value  of  the 

heat  flux  set  to  0.18  K  m  s't,  and  is  the  current  simulated  model  time  in  minutes.  The 
surface  heat  flux  is  then  obtained  by  multiplying  Q„  by  the  surface  air  density  to  give 
units  of  K  kg  m  -  s  *.  The  maximum  downward  heat  flux  at  night  was  set  to  0.005  K  m 
s-i  («  6  W  m-2). 
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Wangara  Day  33  Prescribed  Heat  Flux 


Time  (minutes) 


Figure  3.C.1:  Prescribed  Wangara  Day  33  Surface  Heat  Flux  From  Andre  et  al,  (1978) 


Using  the  prescribed  heat  flux  for  Wangara,  equations  (3.C.2)  and  (3.c.2a)  become 
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3.d.  Conservation  of  Mixing  Ratio 


The  time  scheme  equations  for  the  mixing  ratio  are  developed  in  the  same  manner  as 
the  equations  for  the  potential  temperature.  A  backward  (implicit)  scheme  is  also  used  to 
represent  equation  (2.a.(4).8)  rewritten  similar  to  (3.C.1).  Then,  when  E>0,  the  mixing 
ratio  time  difference  equation  is 
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When  E<0,  0.d.l)  becomes 
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Here,  the  Louis  (1979)  heat-moisture  surface  flux  parameterization  was  used  for  both  the 
Wangara  and  ocean  simulations. 

3.e.  Turbulence  Kinetic  Energy 


The  time  difference  form  of  equation  (2.a.(5).24)  is  written 
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This  is  a  backward  (implicit)  finite  difference  equation  just  as  was  used  for  the  potential 
temperature  and  mixing  ratio  equations,  however,  it  is  also  a  cubic  equation  whose  three 
roots  (one  always  real,  other  two  complex  conjugates  or  real)  are  equal  to  .  When 

the  model  solves  this  cubic  equation  and  the  square  of  the  solution  that  is  always  real 
is  assigned  to 


When  the  entrainment  rate  is  less  than  zero,  the  TKE  at  the  next  time  step  is 
determined  by  applying  a  backward  (implicit)  scheme  to  the  top  equation  of 

C2.a.(6)(b).4),  =  where  D  =  Here  the 

*  (4,)J 

dissipation  is  written  in  partially  implicit  form  so  the  finite  difference  scheme  can  be 
solved  without  using  a  cubic  equation.  As  with  the  forward  scheme  used  for  the 
momentum  equations,  this  partial  implicit  representation  still  provides  a  stable  solution. 
Thus,  the  equation  for  the  TKE  with  E<0  is 
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4.  Simulations 


A  short  description  of  the  Wangara  dataset  is  provided  in  the  next  section.  The 
Wangara  data  is  used  to  validate  the  model  and  study  the  clear  convective  boundary 
layer.  The  last  section  in  this  chapter  briefly  describes  the  ocean  experiment.  The  ocean 
simulation  is  designed  to  study  steady-state  conditions  in  the  PBL. 

4.a.  Land  Simulation 

The  Wangara  dataset  was  compiled  by  Clarke  et.  al.  (1971).  It  consists  of  44  days  of 
boundary  layer  data  from  15  July  to  27  August,  1967.  The  data  was  obtained  from  the 
area  around  Hay,  Australia  located  at  34°30'S,  144°56'W.  The  data  collection  project  was 
given  the  name  "Wangara",  which  means  "west  wind".  Day  33  of  the  Wangara  dataset 
was  used  for  the  land  simulation. 

Day  33  was  characterized  by  clear  skies,  negligible  advection  of  heat  and  moisture, 
and  high  pressure.  The  nearest  front  was  over  1000  km  away.  These  conditions  proved 
perfect  for  study  of  the  clear  convective  boundary  layer.  This  particular  day  has  been 
widely  used  in  boundary  layer  studies  because  of  these  ideal  conditions  and  the  readily 
available  data. 

The  data  includes  temperature  and  mixing  ratio  soundings  every  three  hours  from 
the  surface  to  2000  meters.  Soundings  of  the  u  and  v  components  of  the  wind  are 
provided  every  hour  from  the  surface  to  2000  meters.  The  resolution  of  this  sounding 
data  is  every  50  meters  from  the  surface  to  1000  meters,  and  every  100  meters  for  the 
remainder.  The  ground  temperature  and  the  geostrophic  wind  are  provided  once  an  hour. 
Clarke  provides  additional  data  that  is  not  used  in  this  model. 

4.b.  Ocean  Simulation 


The  data  required  is  minimal  since  steady-state  solutions  are  sought  for  the  ocean 
simulation.  A  constant  sea  surface  temperature  (SST)  is  specified.  The  surface  mixing 
ratio  is  computed  based  on  the  SST.  The  mixing  ratio  at  the  top  of  the  PBL  is  fixed.  The 
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initial  air  temperature  is  specified  to  provide  a  positive  surface  heat  flux  at  the  start  of  the 
simulation.  The  surface  winds  and  the  geostrophic  winds  are  set  to  constants.  The 
potential  temperature  and  winds  at  the  top  of  the  PBL  are  determined  based  on  their 
surface  values  and  constant  lapse  rates.  Finally,  a  divergence  is  specified  to  balance  the 
entrainment  rate  in  the  PBL  depth  prediction  equation. 
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5.  One-Layer  Model  Prognostic  Results 

5.a.  Wangara  Experiments 
5.a.(l)  Twenty-four  Hour  Simulation 

A  24-hour  simulation  using  the  Wangara  Day  33  data  was  run  to  predict  the 
diurnal  variation  of  the  prognostic  variables.  The  dissipation  fraction  (see  Table 
2.a.(6).  1)  and  the  fraction  of  TKE  production  due  to  the  local  rate  of  change  of  TKE 
when  the  entrainment  rate  is  less  than  zero  (Section  2.a.(6)(b))  were  set  to  0.90.  The 
simulation  was  started  at  0900L  with  a  time  step  of  60  seconds.  A  cooling  rate  of  2° 
day^  was  applied  to  the  predicted  mixed  layer  potential  temperature.  The  initial  PEL 
depth  was  set  to  18  mb  (=120  meters).  The  initial  TKE  was  set  to  0.2  m^  s'^.  The 
Coriolis  parameter,  f,  is  equal  to  -8.26  10*^  s  ^  for  Wangara. 

Figure  5.a.(l).l  shows  the  diurnal  change  of  the  PEL  depth.  The  abscissa 
indicates  the  number  of  minutes  into  the  simulation  after  the  start  time. 
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Wangara  PBL  Depth 


Time  (minutes) 

Figure  5.a.(l).l:  Predicted  Diurnal  Az„  for  Wangara  Day  33. 

This  profile  is  typical  of  a  clear  convective  boundary  layer  (CBL).  At  the  start  of  the 
simulation  during  the  early  morning  a  strong  inversion  exists  just  above  the  surface.  The 
boundary  layer  is  shallow  at  this  time  (=100  meters).  The  strong  inversion  present  during 
the  early  morning  acts  to  suppress  the  buoyancy.  Since  buoyancy  is  the  driving  force  in 
CBLs,  the  boundary  layer  grows  slowly  during  this  initial  stage. 

As  the  surface  heating  increases,  the  lapse  rate  transitions  from  stable  to  unstable. 
The  air  just  above  the  surface  warms  enough  to  remove  the  existing  low-level  inversion. 
Figure  5.a.(  1  ).2  shows  d„  -  6s_  and  dg^  -  .  The  first  difference  is  a  measure  of  the 

strength  of  the  surface  inversion  and  the  second  difference  the  strength  of  the  PBL  top 
inversion. 
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Time  (minutes) 

Figure  5.a.(l).2:  Diurnal  d„  -  6^^  and  0^+  -  for  Wangara  Day  33. 

The  surface  heating  becomes  strong  enough  to  remove  the  surface  inversion  after  about 
100  minutes.  This  marks  the  second  stage  when  rapid  boundary  layer  growth  takes  place. 
At  this  time  strong  heating  at  the  surface  creates  buoyant  thermals  which  rise.  The  near¬ 
surface  lapse  rate  is  now  superadiabatic  which  results  in  an  unstable  boundary  layer.  This 
allows  the  thermals  to  continue  to  rise  until  they  reach  the  inversion  marking  the  present 
height  of  the  PBL.  The  large  amount  of  buoyancy  at  this  time  of  day  creates  vigorous 
mixing,  hence  the  name  mixed  layer.  This  causes  the  conservative  variables  to  become 
nearly  uniform  with  height  in  the  mixed  layer.  The  predicted  diurnal  variation  of  the 
mixed  layer  prognostic  variables  are  shown  in  Figures  5.a.(l).3-5.a.(l).6. 
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Wangara  Mixed  Layer  Potential  Temperature 
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Figure  5.a.(l).3:  Predicted  Diurnal  On,  for  Wangara  Day  33. 


Wangara  Mixed  Layer  Mixing  Ratio 
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Figure  5.a.(l).4:  Predicted  Diurnal  qn,  for  Wangara  Day  33. 
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Wangara  Mixed  Layer  Horizontal  Velocity 
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Figure  5.a.(l).5:  Predicted  Diurnal  |v„|  for  Wangara  Day  33. 


Wangara  Mixed  Layer  Turbulence  Kinetic  Energy 
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Figure  5.a.(l).6:  Predicted  Diurnal  em  for  Wangara  Day  33. 

As  the  morning  progresses  and  the  boundary  layer  becomes  deeper,  growth 
occurs  not  only  due  to  buoyancy,  but  also  because  warm  free-atmosphere  air  is  entrained 
into  the  top  of  the  PEL.  This  air  is  mixed  by  the  turbulence  within  the  PEL  causing  the 
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PBL  to  grow.  Entrainment  arises  because  of  penetrative  convection.  This  is  illustrated  in 
Figure  5.a.(l).7. 


e(z) 


Figure  5.a.(l).7:  Illustration  of  the  Process  of  Penetrative  Convection  (Stull,  1991). 

An  air  parcel  in  the  mixed  layer  that  is  initially  warmer  than  the  mean  potential 
temperature  is  positively  buoyant,  and  thus  rises  through  the  layer.  At  this  point  the 
parcel  does  not  require  any  forcing  to  rise.  This  is  free-convection  where  the  parcel  gains 
momentum  during  its  trip  upward.  When  the  parcel  reaches  the  top  of  the  PBL,  it 
encounters  warmer  air  due  to  the  inversion  that  marks  the  transition  from  mixed  layer  to 
free-atmosphere.  The  parcel  then  becomes  negatively  buoyant,  but  continues  to  rise  into 
free-atmosphere  because  of  its  momentum.  This  overshooting  is  called  penetrative 
convection. 

Once  the  air  parcel  has  lost  its  momentum  it  sinks  back  into  the  mixed  layer.  The 
parcel  carries  along  non-turbulent,  warm,  free-atmosphere  air  on  the  return  trip.  The 
positively  buoyant  free-atmosphere  air  becomes  mixed  by  the  turbulence  in  the  mixed 
layer  before  it  has  a  chance  to  escape.  This  capture  and  subsequent  mixing  of  warm  free- 
atmosphere  air  is  the  process  of  entrainment.  Since  less  turbulent  air  is  entrained  into 
more  turbulent  air,  entrainment  only  occurs  in  one  direction  —  down  into  the  PBL. 

Mechanical  mixing  caused  by  wind  shear  at  the  surface  and  top  of  the  PBL  also 
causes  PBL  growth,  but  this  process  is  less  important  in  a  clear  CBL  over  land.  This  is 
shown  by  Figure  5.a.(l).8. 
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Wangara  Diurnal  Shear  and  Buoyancy 


Figure  5.a.(l).8:  Diurnal  B  and  S  Wangara  Day  33. 

During  the  rapid  PBL  growth  period,  buoyancy  production  is  an  order  of  magnitude 
larger  than  shear  production.  Both  buoyancy  and  shear  production  in  this  figure  take  into 
account  the  contribution  due  to  entrainment  when  E>0.  Shear  production  becomes 
important  at  night  when  buoyancy  production  is  negative.  Shear-generated  turbulence 
may  cause  the  nocturnal  boundary  layer  to  grow. 

The  predicted  diurnal  change  in  the  entrainment  rate  is  shown  in  Figure  5.a.(l).9. 
By  mid  morning  when  the  PBL  has  rapidly  grown  to  about  1  km,  the  entrainment  rate  has 
increased  dramatically.  This  gives  an  indication  that  entrainment  is  an  important 
mechanism  for  boundary  layer  growth. 
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Figure  5.a.(l).9:  Predicted  Diurnal  E  Wangara  Day  33. 


The  rapid  decrease  in  PEL  growth  marks  the  third  stage  of  the  diurnal  transition 
of  the  mixed  layer.  At  this  point  rising  thermals  meet  resistance  upon  reaching  the  base 
of  the  inversion  at  the  top  of  the  PEL.  The  inversion  has  increased  in  strength  as 
indicated  in  Figure  5.a.(l).2  which  makes  it  more  difficult  for  penetrative  convection  to 
occur.  Euoyancy  production  is  no  longer  as  effective  in  a  deep  boundary  layer  as  it  was 
when  the  PEL  was  shallow.  The  boundary  layer  continues  to  grow,  however  growth  is 
much  slower.  As  Figure  5.a.(l).9  shows,  the  entrainment  rate  rapidly  drops  off  by  early 
afternoon  which  coincides  with  the  much  slower  growth  rate  of  the  PEL  during  this 
period. 

The  final  stage  in  the  transition  of  the  mixed  layer  occurs  around  sunset.  With  the 
loss  of  daytime  heating,  buoyancy  production  rapidly  approaches  zero.  This  marks  the 
decay  of  turbulence  in  the  mixed  layer.  The  TKE  is  no  longer  maintained  by  buoyant 
production,  and  is  rapidly  dissipated.  The  mixed  layer  then  becomes  decoupled  from  the 
surface.  Since  the  sum  of  the  buoyancy  and  shear  is  now  less  than  the  dissipation,  the 
entrainment  rate  is  allowed  to  become  negative.  This  has  the  effect  of  "crashing"  the 
mixed  layer. 
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The  mixed  layer  depth  decreases  to  its  preset  minimum  of  10  meters  shortly  after 
sunset.  The  mixed  layer  potential  temperature  decreases  continuously  until  sunrise  due  to 
the  constant  downward  heat  flux  and  a  constant  prescribed  radiative  cooling.  The  TKE 
decreases  at  sunset  and  remains  at  its  prescribed  minimum  of  1  10  -  m-  S'2  during  the 
night. 


With  negative  buoyancy  production  and  insufficient  shear  production,  the 
entrainment  rate  remains  negatu  e,  but  it  approaches  zero  after  sunrise  when  the 
buoyancy  production  becomes  positive.  The  boundary  layer  is  expected  to  grow  at  night 
due  to  shear  generated  turbulence  and  other  factors.  There  appears  to  be  a  problem  with 
the  negative  entrainment  parameterization  because  it  does  not  allow  PBL  growth  during 
the  night. 

5.a.(2)  Seventy-two  Hour  Simulation 

The  model  was  then  mn  for  72  hours  to  test  the  response  to  repeat  use  of  the 
Wangara  Day  33  data.  It  was  expected  that  the  prognostic  variable  profiles  would  look 
very  similar  from  day-to-day.  Slight  variations  were  considered  acceptable  because  the 
initial  conditions  at  model  start  time,  0900L  Day  1,  would  not  be  the  same  as  the 
predicted  conditions  24  hours  later,  0900L  Day  2.  These  predicted  conditions  could  be 
considered  the  "new"  initial  conditions  at  the  start  of  the  second  day.  Figure  5.a.(2).l 
shows  the  mixed  layer  PBL  depth  as  a  representative  profile.  The  profile  is  consistent 
from  Day  1  through  Day  3.  Although  not  shown,  the  other  prognostic  variables  were  also 
consistent  throughout  the  simulation. 
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Wangara  PBL  Depth  During  3  Day  Simulation 


Time  (minutes) 

Figure  5.a.(2).l;  Predicted  3  Day  Reusing  Wangara  Day  33  Data  Each  Day. 


5.b.  Ocean  Experiment 


The  ocean  experiment  was  designed  to  obtain  steady-state  solutions  since  no 
database  was  used  for  this  simulation.  The  initialization  of  the  prognostic  variables  is 
detailed  in  Table  2.b.l .  Constants  required  to  initialize  the  prognostic  variables  are  listed 
in  Table  2.b.2.  A  100  hour  simulation  was  run  to  allow  the  variables  to  reach 
equilibrium.  The  PBL  depth  prediction  equation  requires  a  non-zero  divergence  to 
balance  a  positive  entrainment  rate  when  equilibrium  has  been  reached.  The  steady-state 
form  of  this  equation  is 


^m  = 


gE 


(S.b.l) 


A  divergence  of  4  10'^  s'^  was  used  for  this  experiment. 

Figure  S.b.l  shows  the  convergence  of  the  PBL  depth  completely  to  its  steady-state 
value  by  100  hours.  At  equilibrium,  the  local  rate  of  change  terms  in  the  prediction 
equations  are  zero.  As  a  check,  the  steady-state  solution  for  one  of  the  prognostic 
variables  can  be  determined  by  setting  the  local  term  in  the  prediction  equation  to  zero. 
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thus  obtaining  an  equation  for  the  variable  in  terms  of  diagnostic  variables.  The  solution 
to  this  equation  should  equal  the  value  of  the  variable  predicted  by  the  model. 


Ocean  Experiment  PBL  Depth 


Time  (minutes) 

Figure  5.b.l:  Predicted  Az„  Over  Ocean. 


Figures  5.b.2-5.b.5  show  the  progression  of  the  other  prognostic  variables  to  then- 
equilibrium  values.  The  steady-state  equation  for  the  potential  temperature  is 


_  Edg^  +  Vds. 

E  +  V 


(5.b.2) 


If  E»V  then  the  mixed  layer  potential  temperature  will  reach  the  temperature  at  the  PBL 
top  in  equilibrium.  This  shows  that  entrainment  dominates.  If  E«V  then  the  mixed 
layer  potential  tempeiature  will  reach  the  temperature  at  the  PBL  surface  in  equilibrium. 
In  this  case,  the  surface  heating  dominates. 

Using  the  values  of  the  potential  temperature  at  the  top  and  surface  of  the  PBL 
(dg^  =  291.93  K  and  =  288.79  K),  and  the  values  of  E  =  4.38  lO'^  kg  m-2  s'l  and  V  = 
2.36  10-2  kg  m-2  s'l  at  t=6000  minutes,  equation  (5.b.2)  gives0„  =  288.79  K.  This 
compares  almost  exactly  with  6^  =  288.75  K  at  t=6000  minutes  from  Figure  5.b.2. 
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Ocean  Experiment  Mixed  Layer  Potential  Temperature 


Figure  5.b.2:  Predicted  6^  Over  Ocean. 
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Ocean  Experiment  Mixed  Layer  Mixing  Ratio 
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Figure  5.b.3:  Predicted  q„  Over  Ocean. 
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Ocean  Experiment  Mixed  Layer  Horizontal  Velocity 


Time  (minutes) 


Figure  5.b.4;  Predicted  Over  Ocean. 


Ocean  Experiment  Mixed  Layer  Turbulence  Kinetic  Energy 
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Figure  5.b.5:  Predicted  Over  Ocean. 

According  to  Figures  5.b.l-5.b.5,  all  the  prognostic  variables  reach  equilibrium 
about  the  same  time.  Schubert  et  al.  (1979)  developed  a  coupled,  convective-radiative, 
boundary  layer  model  and  performed  several  ocean  simulations  where  they  varied  the  sea 
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surface  temperature  (SST),  the  divergence,  or  both.  In  one  experiment,  the  SST  was 
increased  instantaneously  from  14°C  to  16°C,  and  the  divergence  was  held  constant  at  4 
10'^  s They  found  the  adjustment  time  for  the  PBL  depth  to  reach  steady-state  was 
about  20  times  as  long  as  for  the  other  prognostic  variables.  They  concluded  that  the 
longer  adjustment  time  was  a  general  feature,  at  least  under  some  typical  eastern  ocean 
situations. 


In  Schubert's  study  an  important  dimensionless  quantity  was  introduced  that 
measured  the  relative  importance  of  surface  transfer  and  mixing  across  cloud  top.  This 
quantity  was  adopted  for  the  present  study,  except  that  the  mixing  was  due  to  entrainment 
of  free-atmosphere  air  only  since  no  cloud  effects  were  included.  This  quantity  can  be 
thought  of  as  an  adjustment  ratio  and  has  the  form 


/!  = 


CrV 


DZg  + 


dZg  ' 
dt 


(5.b.3) 


where  Ct  is  the  surface  transfer  coefficient,  V  the  surface  wind  speed,  D  the  divergence, 
and  ZB  the  height  of  the  PBL  in  meters.  If  surface  transfer  dominates  then  the  ratio  is 
large  (about  4  or  5).  The  surface  forcings  rapidly  adjust  the  thermodynamic  variables, 
while  the  slow  mixing  at  the  PBL  top  causes  the  PBL  depth  to  adjust  slowly.  If  the 
mixing  at  the  PBL  top  dominates  then  the  ratio  is  small  (<1).  In  this  case,  the  PBL  depth 
adjusts  in  about  the  same  time  as  the  thermodynamic  variables. 

The  value  of  A  for  the  ocean  experiment  is  shown  in  Figure  5.b.6. 


66 


Ratio  of  Surface  Transfer  to  PEL  Top  Mixing 
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Figure  5.b.6:  A  for  Ocean  Experiment. 

The  ratio  was  never  less  than  about  4.5  which  would  indicate  that  the  PEL  depth  takes 
much  longer  to  adjust  than  the  other  prognostic  variables. 


This  discrepancy  is  resolved  by  comparing  the  entrainment  rate  in  the  ocean 
experiment  with  the  one  used  in  the  Schubert  study.  Figure  5.b.7  depicts  E  for  the  ocean 
experiment. 
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Ocean  Experiment  Entrainment  Rate 


Figure  5.b.7:  E  for  Ocean  Experiment. 


At  equilibrium  the  entrainment  rate  was  small,  but  for  the  initial  portion  of  the  simulation 
the  entrainment  rate  became  very  large.  However,  in  Schubert’s  study  the  entrainment 
rate  remained  at  a  constant  small  value  for  the  entire  simulation.  The  parameterization 
for  E  used  in  the  present  experiment  caused  E  to  become  large  enough  so  that  the  PBL 
depth  adjusted  rapidly. 


If  one  assumes  that  the  equilibrium  value  of  E  obtained  in  the  ocean  experiment  was 

the  value  of  E  for  the  entire  simulation,  then  the  adjustment  time  for  the  PBL  is  obtained 

HArt 

by  solving  the  differential  equation  for  the  PBL  depth,  •  The 

solution  to  this  equation  with  the  divergence  and  E  constant  is 


V  •  V  ^ 


.+(4Pm)o 


(5.b.4) 


where  is  the  e-folding  time  (time  for  variable  to  decrease  to  1/e  of  its  original 
value)  and  is  the  initial  PBL  depth.  As  oo  the  PBL  depth  reaches  its 

equilibrium  value  of  Ap„  Equation  (5.b.4)  can  be  manipulated  to  get  a  relation 


for  t 


€•-  fold  * 
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This  relation  is 


M  4p„- 


gE 

V  •V 


m  J 


(4>Jo- 


gE 


^  enfold 


-V«  V. 


(5.b.5) 


The  e-folding  time  obtained  with  E  at  its  steady  state  value  of  4.382529  lO'^  kg  m-2  s-i,  a 

divergence  of  4  10’^  s  ^  an  initial  depth  of  5817.8  Pa  and  a  final  depth  of  10736.5  Pa  was 
25.6  days.  The  adjustment  time  is  approximately  3  times  which  is  about  77  days. 

This  is  about  18  times  as  long  as  the  PBL  depth  actually  took  to  adjust  in  the  ocean 
experiment  which  corresponds  excellently  with  the  Schubert  study.  This  shows  that  the 
differences  in  E  between  the  present  study  and  Schubert’s  study  are  the  key  to  the  rapid 
adjustment  of  the  PBL  in  the  present  study. 
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6.  One-Layer  Model  Diagnostic  Discussion 

This  chapter  provides  a  brief  overview  of  diagnostic  variables  determined  in  the 
model  following  Randall  et  al.  (1992).  These  variables  include  the  fractional  area 
covered  by  rising  motion,  a,  the  convective  mass  flux.  Me,  plume-scale  variance 
transport,  pw'y/'y/'  (where  \|/  is  an  arbritary  scalar  such  as  the  potential  temperature  or 
water  vapor  mixing  ratio),  value  of  \\f  at  levels  S  and  B  for  upward  and  downward  moving 
parcels,  (V|/u  or  d)s  or  dissipation  time  scale,  Xdis,  dissipation  rate  of  \|/  at  levels  S  and  B, 

(Vdis)s  or  Vertical  gradient  of  surface  transfer  coefficients,  Q  and  Q., 

az 

Richardson  number,  and  Richardson  number  limits. 

6.a.  Convective  Mass  Flux  Model 

The  scalar,  \|/  satisfies  the  conservation  equation 

^i^  =  -V«(pvv^)-^(pwvr)  +  S  ,  (6.a.l) 

at  az 

where  the  local  change  and  the  del  operator  are  defined  on  constant  height  surfaces,  and 
Sv|r  is  the  source  of  \|/  per  unit  mass  per  unit  time.  The  area  average  of  the  scalar  is  given 
by 

¥  =  +  (6.a.2) 

and  the  upward  turbulent  flux  of  \|/  is 

=  pVijP  =  p[{w^  -  vv)(v'„  -  W)(^  +  -¥)U-  cy)] 

or  (6.a.3) 

where  the  convective  mass  flux  is  defined  as 
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Me  =pa{l-G){w^-Wj). 


(6.a.4) 


The  convective  mass  flux  can  not  be  determined  by  the  model  using  (6.a.4)  because  the 
vertical  velocities  of  upward  and  downward  moving  parcels  are  not  known  nor  predicted. 
The  convective  mass  flux  can  also  be  written  in  terms  of  the  fractional  area  covered  by 
rising  motion  and  the  turbulence  kinetic  energy.  The  former  is  diagnostically  determined 
(see  Section  6.c)  using  the  entrainment  rate  and  ventilation  mass  flux  which  are 
calculated  by  the  model,  and  the  latter  is  predicted  by  the  model.  The  definition  (6.a.4)  is 
useful,  however,  in  developing  an  equation  for  the  plume-scale  variance. 

The  plume-scale  variance  transport  is 

pw'y/Y  =  p[cr(w„  -  vv)(v^„  -  +  (7  -  -  w){\i/,  -  ^)'] , 

or  using  (6.a.2)  where  w=\|/. 


pw'y/Y  =  pa{l-  a){l  -  2a){w^  -  w,){y/^  -  y/,)\ 


or  using  (6.a.3)  and  (6.a.4), 


pw'Yw'  =  {l-2o) 


(M. 

Me 


(6.a.5) 


Equations  for  v|/u  and  \|/d  are  obtained  by  substituting  (6.a.3)  into  (6.a.2)  after 
rearrangement  which  gives 


and 


(6.a.6) 


(6.a.7) 
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6.b.  Matching  Convective  Mass  Flux  with  Ventilation  and  Entrainment  Mass  Flux 
For  the  ventilation  layer,  the  bulk  aerodynamic  formula  used  by  the  model  is 

=  (6b.l) 

The  ventilation  mass  flux  can  be  matched  to  the  convective  mass  flux  at  the  top  of  the 
ventilation  layer  (at  level  S)  with  the  following  assumptions:  (1)  The  fluxes  at  the  top  of 
the  ventilation  layer  are  entirely  due  to  convective  circulations,  and  the  small-eddy  fluxes 
are  negligible  at  S.  This  is  a  typical  assumption  in  the  boundary  layer  where  the  small 

eddies  are  important  very  near  the  surface  (viscous  dissipation.  Re  =  =  1,  where  U  is 

V 

the  horizontal  velocity,  L  is  the  length  scale  of  the  eddy,  and  v  is  the  viscosity),  but  in 
most  of  the  surface  layer  the  Reynolds  number  is  large  (since  U  and  L  are  large  and  v  is 
small  compared  to  their  values  in  the  viscous  sublayer)  and  viscous  effects  are  no  longer 
important.  (2)  The  ventilation  layer  is  thin  (the  model  assumes  the  ventilation  ana 
entrainment  layers  are  infinitesimal).  The  ventilation  mass  flux  can  then  be  matched  to 
the  convective  mass  flux  at  level  S  with  these  assumptions  and  equation  (6.a.3).  This 
gives 


Hvs-  -¥s)  =  McAVu  -  ¥,)s-  (6-b.2) 

Since  the  small  eddies  are  important  near  the  surface,  they  will  dilute  air  that  rises 
from  the  surface  and  air  that  descends  from  the  interior  of  the  PBL.  To  account  for  this 
mixing,  a  mixing  parameter,  Xv»  is  used  so  that 

{¥u)s  -¥s=  XviVs-  ~Vs)^  (6.b.3) 

where  0<Xv-^-  When  the  mixing  parameter  equals  1,  no  mixing  occurs  by  the  small 
eddies  and  =  y^s--  Mixing  by  the  small  eddies  increases  as  the  mixing  parameter 
decreases  from  1.  Using  (6.a.2),  (6.b.2)  and  (6.b.3)  results  in 

Mr  c 

=  (6.b.4) 
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A  similar  matching  of  the  fluxes  in  the  entrainment  layer  at  level  B  leads  to 


Xe- 


M. 


C.B 


=  o 


8- 


(6.b.5) 


Now,  if  Me  and  a  are  assumed  to  be  independent  of  height  then  (6.b.4)  and  (6.b.5)  can  be 
combined  to  give 


(7  = 


1  +  -^' 


y  Xe 


(6.b.6) 


An  equation  for  Me  in  terms  of  the  entrainment  rate,  the  ventilation  mass  flux,  and  the 
mixing  parameters  is  obtained  by  inserting  (6.b.6)  into  (6.b.4)  or  (6.b.5)  which  gives 


"  (EIXrMVIXvY 


(6.b.7) 


The  model  does  not  determine  a  or  Me  using  equations  (6.b.6)  and  (6.b.7)  because  it 
does  not  contain  a  parameterization  for  the  mixing  parameters.  The  next  section  presents 
an  equation  for  Me  in  terms  of  a  and  the  TKE.  This  equation  is  equated  with  (6.b.7)  to 
deduce  a  parameterization  for  the  mixing  parameters  where  they  are  equal  to  the  same 
quantity.  The  parameterization  is  not  applied  directly  by  the  model,  but  is  used  to 
simplify  (6.b.6). 


6.C.  Diagnostic  Equations  for  Me  and  o  Using  the  TKE 


Assuming  the  density  of  air  is  approAimately  constant  with  height  in  the  PBL  (since 
the  PBL  depth  is  typically  only  1-2  km),  the  vertically  averaged  TKE  (e„i)  is  related  to 
the  variance  of  the  vertical  velocity  by 


(6.C.1) 


73 


where  33=0.316.  This  equation  is  simply 


7 

Kinetic  Energy  _  2 


area 


area 


Now,  the  variance  of  the  vertical  velocity  is  written 


vv'-  =tj(/-a)(w„ 


or  using  (6.a.4), 


H’  = 


Me 


(6.C.2) 


All  the  quantities  on  the  right  hand  side  of  (6.C.2)  are  assumed  to  be  independent  of 
height.  Then,  substituting  this  equation  into  (6.c.l)  and  integrating  gives 


o,e„,p„,Az  = 


/  Ml 

2  pia{I-o) 


Pm^  ’ 


or 


(6.C.3) 


Once  the  final  equation  for  a  is  determined  then  Me  can  be  calculated  using  (6.C.3). 


Setting  equation  (6.C.3)  and  (6.b.7)  equal  to  each  other  results  in 


(6.C.4) 


Then  substituting  for  a  using  (6.b.6)  to  obtain 


EV 


2p:.a,e„, 


(6.C.5) 


A  plausible  parameterization  for  the  mixing  parameters  based  on  (6.C.5)  is  then 
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Xv  =  Xe 


'  2p^  a,e 

r  m^3  n 


(6.C.6) 


Finally  with  this  parameterization,  (6.b.6)  reduces  to 


V 


(6.C.7) 


The  model  calculates  the  entrainment  rate  and  the  ventilation  mass  flux,  and  then  (6.C.7) 
and  (6.C.3)  are  used  to  determine  a  and  Me- 

6.d.  PEL  Interior  Diagnostics 


The  balance  for  the  variance  of  \\f  in  the  PEL  interior  is  written 


3  fZ  3\j/  I  S  !  /  /  ^ 


(6.d.l) 


where  the  local  change  of  the  variance  is  due  to  production  of  variance,  vertical  ti.'  spos 
of  variance,  and  dissipation  of  variance  (see  (6.d.3)).  Advection  by  the  mean  flow  has 
been  ignored  and  \|/  is  assumed  to  be  a  conservative  variable.  The  variance  is  given  by 

or  using  (6.a.3), 


Mr 


(6.d.2) 


The  triple  correlation  portion  of  the  triple  correlation  term  is  just  the  plume-scale  variance 
transport  (6. a. 5).  The  dissipation  rate  used  by  the  model  is 


y/''  _  (j{I-o)(  /y 
[  Mr 
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(6.d.3) 


o{l  -  o)(l  -  2oy  \ 

(  F 

f  ' 

[Mc 

where 


^dis  ~ 


il-2ay 


(6.d.4) 


Equation  (6.d.4)  is  used  to  calculate  the  dissipation  time  scale  for  tj/  based  on  a  and  the 
parameter,  f ,  which  is  set  during  a  model  simulation.  The  model  determines  the 
dissipation  rate  of  the  variance  of  \|/  using  the  lower  equation  of  (6.d.3). 


The  last  diagnostic  to  be  determined  in  this  section  is  the  vertical  gradient  of  V|/. 
Writing  (6.d.l)  using  (6.d.2),  the  plume-scale  variance  transport,  and  the  top  equation  of 
(6.d.3)  gives 


dt 


f  F  ^ 

0(1-0) 

J 


p  ^z 


Mc(l-2o)  I  d 
o(I-o)  pdz 


(  F^ 

o(l-o) 

V 

2o(l-o) 


''dis 


F 

V 


(6.d.5) 


An  equilibrium  solution  to  (6.d.5)  can  be  found  by  setting  the  local  time  derivative  to 
zero.  The  equation  then  contains  a  first  order  derivative  in  z  which  requires  only  a  single 
boundary  condition  to  solve.  The  boundary  condition  is  applied  at  level  S  if  a<l/2 
(boundary  layer  driven  by  surface  heating),  and  at  level  B  if  a>l/2  (boundary  layer 
driven  by  entrainment).  To  satisfy  both  (6.b.l)  (surface  flux)  and 

an  additional  condition  must  be  specified.  Choosing  to  be  constant  with  height  will 

az 


force  the  differential  equation  to  be  satisfied  at  both  boundaries. 


Using  the  hydrostatic  equation  and  the  conditions  above,  (6.d.5)  becomes 


{I-2o)(  dF^ 
M(.  (  dty 


F^  _dy/ 

SfJ. ,  diJ  ' 


(6.d.6) 
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(6.d.7) 


where  dp.  = 


gMci 


a{l-o){l-2a) 


Then,  using  the  surface  and  top  fluxes  as  boundary  conditions,  the  solution  of  (6.d.6)  is 


1-exp 


P  Pb 


8p. 


exp\ 


f  \ 

P-Pb 


dp. 


-exp 


J 


1  -  exp\ 


dp. 


V  yj 


(6.d.8) 


where 


Sy/  _  {l-2a) 


dp  M(~dp. 


dp. 


1-exp 


^dp.  j 


(6.d.9) 


Equation  (6.d.9)  is  used  by  the  one-layer  model  to  obtain  the  vertical  gradient  of  \|/. 

Assuming  a  is  close  to  1/2  and  using  the  binomial  expansion,  equations  (6.d.8)  and 
((6.d.9))  become 


Ps-P 


('v), 


P-Pb 


4p^ 


and 


P-Pb 

\ 


8t: 


f . 


Ps-P 

y 


[('vl-('v)J' 


(6.d.l0) 


dy/  ^  _ _ 

dp  G{l-a){ApJ 


A  dp. 


Equation  (b.d.lO)  is  an  approximation  to  (6.d.8)  keeping  first  order  terms.  Equation 
((6.d.l  1)  is  an  approximation  to  (6.d.9)  keeping  second  order  terms.  T'^en,  if 

» |(^v')fi|’  '''hich  is  typically  true  in  a  convective  boundary  layer,  ((6.d.l  1)  reduces 

to 
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The  vertical  gradient  of  Q  is  also  determined  by  (6.d.l2)  in  the  one-layer  model. 

6.e.  Surface  Transfer  Coefficient  Using  TKE 

The  bulk  aerodynamic  formula,  -y/s)^  can  be  written  by  specifying 

the  ventilation  mass  flux  (V)  using  the  surface  density,  surface  wind  speed,  and  a  surface 
transfer  coefficient 


^  =  PsCtK[  (6-e.l) 

Based  on  Randall  and  Shao  (1990),  the  ventilation  mass  flux  can  be  related  to  the  TKE  by 
^  =  PsCr4^-  (6-e.2) 


In  (6.e.2),  the  square  root  of  the  TKE  is  'acting'  as  the  velocity.  Since  turbulent  flux 
requires  TKE  and  V  is  a  measure  of  this  flux  at  the  surface,  it  seems  reasonable  that  V  is 
propoitional  to  Another  reason  to  favor  (6.e.2)  over  (6.e.  1)  is  that  turbulence  can 

occur  in  the  absence  of  a  mean  wind  (i.e.,  when  there  is  positive  buoyancy  production). 
As  long  as  TKE  exists,  (6.e.2)  will  determine  V  regardless  of  the  value  of  the  mean  wind 
Both  Ct  and  Q.  are  determined  by  the  model  using  (6.e.l)  and  (6.e.2)  respectively. 


6.f.  Richardson  Number  and  Limits 


The  Richardson  number  is  determined  using  the  equation  listed  in  Section  2.a.(6)(a). 
This  equation  is 


Kl  — - 7— - 


(e.). 


(6.f.l) 


When  the  inversion  is  strong,  Ri»l  then 


(6.f.2) 


where  Cp  is  the  specific  heat  of  air  at  constant  pressure. 


When  there  is  no  inversion,  Ri=0  (see  equation  (2.a.(6)(a).7))  then 


Urn  = 

Ri=0 


E 


=  1. 


(6.f.3) 


The  strong  and  no  inversion  limits  are  determined  by  the  model  using  (6.f.2)  and  (6.f.3) 
respectively. 


79 


7.  One-Layer  Model  Diagnostic  Results 

Results  are  provided  covering  the  convective  growth  period  of  the  Wangara  Day  33 
simulation.  This  period  is  roughly  from  0900L  to  1600L,  and  includes  rapid  growth  of 
the  PEL  during  the  mid-morning  and  slower  growth  during  the  afternoon.  The  results 
point  out  the  importance  of  buoyancy  and  entrainment  in  the  growth  of  a  clear  convective 
PEL  when  the  PEL  top  is  below  a  weak  or  non-existent  inversion.  During  the  afternoon 
when  the  inversion  is  strong,  surface  heating  is  still  significant  which  continues  creating  a 
large  amount  of  buoyancy,  but  this  buoyancy  is  largely  ineffective  in  penetrating  the 
inversion  layer.  The  strong  inversion  layer  also  limits  the  entrainment  rate.  The  small 
amount  of  entrainment  present  is  largely  balanced  against  subsidence,  hence  the  PEL  is 
quasi-steady-state  during  the  afternoon.  Results  are  also  shown  for  the  steady-state  ocean 
experiment. 

7.a.  Wangara  Results  for  the  Fractional  Area  Covered  by  Rising  Motion 

Figure  7.a.l  shows  the  fractional  area  covered  by  rising  motion,  a,  as  a  function  of 

time. 
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Fractional  Area  Covered  by  Rising  Motion 


Figure  7.a.  1 :  a  for  Wangara  Day  33  0900-  1600L. 

By  mid-morning,  a«l  which  is  when  rapid  PBL  growth  is  occurring.  After  1200L,  a 
increases  steadily  as  convective  growth  begins  to  diminish.  The  fractional  area  exceeds 
0.5  after  1520L.  At  this  time  convection  is  no  longer  significantly  affecting  the  PBL 
depth. 

In  Figure  7.a.2  the  plume-scale  variance  transport  of  the  potential  temperature  at 
levels  S  and  B  has  been  overlaid  with  o. 
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Overlay  of  Variance  Transports  and  Fractional  Area 


H 
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Figure  7.a.2:  and  c  for  Wangara  Day  33  0900-1600L. 


It  is  clear  from  the  figure  that  the  plume-scale  variance  transport  of  0  at  level  S  dominates 
during  rapid  convective  growth  when  a«l.  While  a<l/2,  >  0  which 

indicates  the  surface  is  transporting  variance  upwards.  When  a  equals  1/2,  both 
(pw'0'0')^  and  (pw'0'0')^  are  zero.  Finally,  when  a  exceeds  1/2,  [pw'd'O'^^ 

and(pw'0'0')^  are  less  than  zero.  At  this  point  the  entrainment  layer  is  exporting 

variance  downward  into  the  PBL.  This  variance  export  balances  subsidence  keeping  the 
PBL  in  a  quasi-steady-state. 

The  convective  mass  flux  is  shown  in  Figure  7. a.3.  The  minimum  occurs  when 
a«l,  and  the  maximum  occurs  when  a=l/2  while  the  TKE  is  still  large.  Me  is  small 
when  the  convection  is  intense  because  M^-  o^-yJa{l  -o).  As  a  increases  and  the  TKE 

decreases  during  the  late  afternoon.  Me  decreases. 


82 


Convective  Mass  Flux 


Time  (minutes) 

Figure  7.a.3;  Me  for  Wangara  Day  33  0900-1600L. 

The  updraft  (u)  and  downdraft  (d)  properties  of  0  and  q  at  level  S  are  depicted  in 
Figures  7.a.4  and  7.a.5. 


Updraft  and  Downdraft  Potential  Temperatures  at  Level  S 
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Figure  7.a.4:  and  (0,)  for  Wangara  Day  33  0900-1600L. 


updraft  and  Downdraft  Mixing  Ratios  at  Level  S 


Time  (minutes) 

Figure  7.a.5:  (<7„)^  and  (<7^)^  for  Wangara  Day  33  0900-1600L. 

These  figures  indicate  that  the  updrafts  are  warmer  and  wetter  than  the  downdrafts.  The 
boundary  layer  is  being  heated  from  the  surface,  and  the  highest  amount  of  moisture  is 
near  the  surface.  Hence,  the  updrafts  which  are  coming  from  a  region  that  is  warm  and 
moist,  should  be  warm  and  wet  compared  to  the  downdrafts  which  come  from  a  relatively 
dry  and  cool  region. 

Initially  the  surface  heating  rate  is  greater  than  the  surface  heat  transport.  Thus, 
rising  air  near  the  surface  heats  rapidly  before  ascending.  This  causes  the  updraft 
potential  temperature  to  increase  rapidly.  Eventually,  the  surface  heat  transport  exceeds 
the  surface  heating.  Also,  the  intense  heating  and  convection  have  removed  some  low- 
level  available  moisture.  The  surface  air  then  rises  before  it  can  be  heated,  and  it  rises  in 
a  region  of  less  moisture.  This  causes  the  updraft  potential  temperature  to  decrease  for  a 
short  period.  Finally,  when  the  convection  becomes  less  intense,  the  heating  rate  again 
exceeds  heat  transport.  The  moisture  loss  also  decreases.  At  this  point  the  updraft 
potential  temperature  begins  to  increase,  but  not  as  rapidly  because  of  less  intense  surface 
heating. 
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The  downdraft  potential  temperature  increases  rapidly  in  the  morning  when  the 
heating  is  intense  and  the  heat  transport  is  rapid.  Heat  is  brought  quickly  into  the  source 
region  of  the  downdrafts.  Initially,  moisture  is  also  brought  into  this  source  region.  In 
the  afternoon,  as  the  surface  heating  decreases  and  moisture  is  carried  away  from  the 
source  region,  the  downdraft  potential  temperature  increases  much  more  slowly. 

Figures  7.a.6  and  7.a.7  show  the  updraft  and  downdraft  properties  at  level  B. 


Updraft  and  Downdraft  Potential  Temperatures  at  Level  B 


Figure  7.a.6:  and  {Oj)^  for  Wangara  Day  33  0900-1600L. 


85 


W) 

6a 


o 

■a 

c3 

0^ 


6a 
c 
•  ^ 

X 
•  ^ 


Updraft  and  Downdraft  Mixing  Ratios  at  Level  B 
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Figure  7.a.7:  {q^)^  and  (<?^)^  for  Wangara  Day  33  0900-1600L. 


These  properties  are  largely  controlled  by  the  inversion  at  the  top  of  the  PBL.  In  the 
morning  and  afternoon  when  the  inversion  is  strong,  updraft  air  at  B  is  cooler  than 
downdraft  air  at  B.  During  the  convective  period  when  there  is  no  inversion,  the  air  from 
below  is  rapidly  heated.  The  updraft  air  at  B  then  comes  from  a  warmer  source  than  the 
downdraft  air  at  B.  The  updraft  air  at  B  is  always  wetter  than  the  downdraft  air.  The 
updraft  mixing  ratio  increases  rapidly  to  a  high  value  for  a  short  time  when  the  vertical 
moisture  transport  is  large  during  convection.  Mixing  brings  this  large  value  back  down. 
Both  the  updraft  and  downdraft  mixing  ratios  decrease  in  the  afternoon  because  the 
sources  of  moisture  from  above  and  below  decrease  due  to  heating  and  mixing. 

7.b.  Wangara  PBL  Interior  Results 

Interior  results  were  obtained  for  the  convective  period  of  Wangara  Day  33  using 
four  different  values  of  f  .  Dissipation  rates  for  0  and  q  and  the  dissipation  time  scale 
were  determined  with  f  set  to  1  second.  These  diagnostics  are  just  f  times  their  values 
at  f  =  /  second  for  other  settings  of  f.  Figure  7.b.l  shows  the  dissipation  time  scale. 
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Figure  7.b.  1 :  for  Wangara  Day  33  0900-1600L. 

The  minimum  in  occurs  during  the  maximum  convection  around  1  lOOL  when  o«l. 
This  is  when  the  surface  heating  is  the  most  intense  and  when  the  smaller  eddies  would 
be  the  most  effective.  As  a  1/2  during  the  afternoon  the  PBL  becomes  more  mixed. 
The  variance  transports  decrease  and  the  time  scale  for  dissipation  increases.  When 
a=l/2  at  1520L  <x>,  hence  the  sharp  peak  in  the  figure. 

Figures  7. b.2  and  7.b.3  contain  and  at  levels  S  and  B. 
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Figure  7.b.2;  (e^)^  and  {ee)g  for  Wangara  Day  33  0900-1600L. 
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Figure  7.b.3:  and  for  Wangara  Day  33  0900-1600L. 

The  dissipation  rates  are  highest  when  the  fluxes  are  the  largest  during  mid-morning  rapid 
growth.  For  the  potential  temperature  the  surface  flux  dominates  over  the  flux  at  B  due  to 
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surface  heating,  thus  the  potential  temperature  dissipation  rate  at  S  is  much  greater  than  at 
B.  For  the  mixing  ratio  the  opposite  is  true.  The  large  entrainment  rate  present  when 
rapid  growth  is  occurring  causes  the  mixing  ratio  flux  at  B  to  be  much  greater  than  at  S. 
All  the  dissipation  rates  approach  zero  as  cr  7/2. 

The  next  set  of  figures  shows  the  gradients  of  6  and  q  with  height  using  f  equal  to 
10,  100,  and  1000  seconds.  The  gradients  for  f  =  7  second  are  not  included  because  they 
are  too  large.  The  gradient  profiles  of  d  are  in  Figure  7.b.4. 


Potential  Temperature  Gradient 


Time  (minutes) 

Figure  7.b.4:  —  with  f  =  10, 100,  and  1000  Seconds  for  Wangara  Day  33  0930-1600L. 
oz 


The  gradient  of  6  with  f  =  10  seconds  seems  reasonable  between  210  and  420  minutes 
based  on  the  actual  Wangara  temperature  profile.  The  other  gradients  look  plausible 
during  the  entire  period,  but  the  gradient  with  f  =  1000  seconds  is  the  most 
representative.  This  is  particularly  true  during  the  convective  growth  period  when  this 
gradient  indicates  the  potential  temperature  is  increasing  with  height.  Observations  have 
verified  that  the  upward  heat  flux  is  countergradient  (Wyngaard  and  Brost  1984).  Based 
on  the  Wangara  data  and  the  gradient  profiles  shown,  f  should  be  between  100  and  1000 
seconds  for  typical  convective  boundary  layers. 

Figure  7.b.5  gives  a  similar  set  of  gradient  profiles  for  the  mixing  ratio. 
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Mixing  Ratio  Gradient 


Time  (minutes) 

Figure  7.b.5  with  f  =  10,  100,  and  1000  Seconds  for  Wangara  Day  33  0930-1600L. 
dz 

Here  again,  the  10  second  profile  is  only  reasonable  during  a  portion  of  the  period.  The 
other  profiles  produce  good  results  all  the  time.  It  would  seem  that  a  f  between  100  and 
1000  seconds  would  work  for  q  as  well.  The  q  gradient  profiles  are  also  consistent  with 
observations  showing  the  mixing  ratio  decreasing  with  height  in  a  convective  PBL. 

The  gradient  of  B  was  also  determined  using  equation  (6.d.l2).  The  gradient  using 
(6.d.l2)  is  independent  of  f  because  the  f  in  the  numerator  of  (6.d.l2)  cancels  out  with 
the  f  in  the  denominator  (part  of  the  dp^  term).  Figure  x  shows  this  gradient. 
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Potential  Temperature  Gradient  Using  Equation  (6.d.l2) 


Time  (minutes) 

Figure  7.b.6:  —  Using  (6.d.  12)  for  Wangara  Day  33  0930-1600L. 
az 

This  profile  looks  reasonable  at  all  times  and  it  shows  the  large  gradient  during  the 
morning  before  the  PEL  has  become  mixed,  rapid  decrease  in  the  gradient  during  the 
mid-morning  convective  period,  and  the  near  zero  gradient  in  the  afternoon  after  mixing 
has  occurred. 

7.C.  Wangara  Surface  Transfer  Coefficients 

Figure  7.C.1  is  a  comparison  of  the  surface  transfer  coefficient  computed  by  using 
the  surface  velocity  with  the  coefficient  calculated  using  the  square  root  of  the  TKE. 
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Surface  Transfer  Coefficients 


Time  (minutes) 


Figure  7.c.  1:  Q  and  Cj,  for  Wangara  Day  33  0900-1600L. 

The  transfer  coefficient  computed  using  is  about  an  order  of  magnitude  larger  than 
the  coefficient  computed  using  |v„|  since  the  square  root  of  the  TKE  is  about  1/10  of  the 

surface  velocity.  The  minimum  occurs  in  this  coefficient  when  em  is  at  its  maximum 
value  from  mid-morning  through  early  afternoon. 


Figure  7.c.2  shows  a  scatter  plot  of  Q.  versus  the  negative  of  a  bulk  Richardson 
number  defined  by 


^huik  — 


(7.C.1) 
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Figure  7.C.2:  Scatter  Plot  of  Q-  Versus  -Rieuik  for  Wangara  Day  33  0900-1600L. 

There  appears  to  be  a  relationship  between  C^.-  and  the  bulk  Richardson  number.  The 
figure  indicates  that  there  are  two  families  of  curves  which  likely  means  that  Cj..  also 
depends  on  another  variable. 

7,d.  Wangara  Calculation  of  Richardson  Number  and  Limits 
A  plot  of  the  Richardson  number  is  shown  in  Figure  7.d.l. 
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Richardson  Number 


Time  (minutes) 


Figure  7.d.l:  Ri  for  Wangara  Day  33  0900-1600L. 


The  Richardson  number  is  zero  during  the  unstable  convective  growth  period  when  there 
is  no  inversion.  At  this  time,  the  limit  when  Ri=0  should  be  1.  As  the  PBL  becomes  well 
mixed  during  the  late  afternoon  the  inversion  strengthens.  The  Richardson  number 
increases  as  a  result.  The  limit  for  Ri»l  should  approach  0.2  by  late  afternoon.  Figures 
7.d.2  and  7.d.3  are  plots  of  the  limits. 
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Richardson  Number  Equal  to  Zero  Limit 
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Figure  7.d.2:  Ri=0  Limit  for  Wangara  Day  33  0900-1600L. 


Richardson  Number  Large  Limit 
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Figure  7.d.3:  Ri»l  Limit  for  Wangara  Day  33  0900-1600L. 

As  indicated  in  the  figure,  the  no  inversion  limit  is  almost  exactly  1  during  the  rapid 
growth  period.  During  the  late  afternoon,  the  Ri»l  limit  does  approach  0.2,  but  it  is  a 
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little  too  small.  This  may  mean  that  there  is  not  an  exact  bal"  e  between  buoyant 
production  and  dissipation  of  TKE  as  assur'^d  in  the  entrainment  closure. 

7.e.  Ocean  Experimc  Fractional  Area  Covered  by  Rising  Motion  Results 

Figure  7.e.l  is  a  for  the  ocean  experiment.  The  initial  difference  between  the  SST 
and  air  temperature  creates  an  upward  surface  temperature  flux.  As  a  result,  <t<1/2  for  a 
short  time.  In  equilibrium,  a  negative  surface  heat  flux  is  required  to  balance  a  positive 
entrainment  rate.  Thus,  in  steady-state,  entrainment  dominates  and  a>l/2.  The  boundary 
layer  would  be  characterized  by  wide  updrafts  with  zones  of  narrow  downdrafts. 


Ocean  Experiment  Fractional  Area  Covered  by  Rising  Motion 


Time  (minutes) 

Figure  7.e.l:  a  for  Ocean  Experiment. 

The  plume-scale  variance  transport  of  the  potential  temperature  at  levels  S  and  B  is 
shown  in  figure  7.e.2. 
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Ocean  Experiment  Potential  Temperature  Variance  Transports 


Time  (minutes) 

Figure  7.e.2:  {^pw'O'd'^^  and  {^pw'd'9'^^  for  Ocean  Experiment. 

When  a<l/2  both  and  {^pw'9'9')^  are  greater  than  zero.  They  transition  from 

negative  to  positive  and  back  to  negative  when  o  becomes  less  than  1/2  and  then  greater 
than  1/2.  Unlike  Wangara,  {^pw'9'9'^^  never  substantially  dominates  over  {^pw'9'9'^^. 

In  steady-state,  the  magnitude  of  (pw'0'0')^  is  greater  than  the  magnitude  of  (^pw'9'9')^. 

Since  the  transports  are  negative,  the  entrainment  layer  exports  variance  into  the  PBL 
which  balances  with  the  dissipation  at  the  surface. 

Figure  7.e.3  shows  the  convective  mass  flux. 
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Figure  7.e.3:  Me  for  Ocean  Experiment. 

The  convective  mass  flux  peaks  when  o<l/2  and  the  TKE  is  large.  This  marks  the  short 
convective  period  when  the  PBL  grows  The  minimum  of  Me  occurs  when  a>l/2  and  the 
TKE  is  at  its  lowest  value.  Here  the  surface  heat  flux  is  negative  and  the  entrainment  rate 
is  at  its  minimum.  For  a  brief  period,  the  divergence  is  removing  mass  faster  than  it  can 
be  replaced  by  entrainment.  There  is  no  convection  with  the  negative  heat  flux  to  aid  in 
PBL  growth.  As  a  result,  the  PBL  depth  levels  off  and  then  decreases  until  'he 
entrainment  rate  increases  sufficiently  to  balance  the  divergence  At  steady-state  the  TKE 
and  Me  are  about  twice  their  minimum  value.s. 


For  Wangara,  Me  was  at  its  minimum  value  during  the  most  intense  convection. 
The  entrainment  rate  was  about  20  times  as  large  as  the  ventilation  mass  flux.  When 
E»V,  equation  6.C.7  can  be  approximated  by 


The  ratio  of  V  to  E,  and  a  become  .small  when  E»V.  This  will  cause  Me  to  be  small 
even  though  vigorous  convection  is  taking  place  and  the  PBL  is  growing  rapidly.  For  the 
ocean  experiment  E  was  only  about  1.2  times  V  during  convection.  The  value  of  a  was 
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less  than  1/2,  but  much  larger  than  the  minimum  from  the  Wangara  simulation.  With  a 
near  1/2  and  the  TKE  large.  Me  was  at  its  maximum  at  the  same  time  as  the  convection. 
Caution  must  be  used  when  comparing  the  convective  mass  flux  to  PBL  growth.  Growth 
may  occur  with  a  low  value  if  the  entrainment  rate  and  ventilation  mass  flux  are  large 
enough  to  balance  subsidence  and  divergence. 

Figures  7.e.4  and  7.e.5  present  the  updraft  and  downdraft  properties  for  the  potential 
temperature  and  mixing  ratio  at  level  S. 


Ocean  Exp)eriment  Up/Down  Potential  Temperatures  at  Level  S 

289  ^ -  -  - : 

i 

288.8  ! 

g  "/  ! 

§  288.6  -f  \  , . j 

I  [:  v;:-— X - - ^ 

I  288.4  #  / 

288.2  'r 


2gg  - 1 - ' - ' - - - ^ - 

0  1000  2000  3000  4000  5000  6000 

Time  (minutes) 


Figure  7.e.4:  and  for  Ocean  Experiment. 
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Ocean  Experiment  Up/Down  Mixing  Ratios  at  Level  S 


Time  (minutes) 

Figure  7.e.5:  {q^)^  and  for  Ocean  Experiment. 

The  updrafts  are  initially  warmer  and  wetter.  The  surface  heat  flux  transports  heat 
vertically  which  warms  the  downdrafts.  Eventually,  the  downdrafts  exceed  the 
temperature  of  the  updrafts.  When  the  heat  flux  becomes  negative,  6^5  begins  to  decrease. 

This  causes  the  downdraft  temperature  to  decrease  despite  the  smaller  positive 
contribution  from  the  negative  heat  flux  (see  equation  (6.d.9)).  The  updraft  potential 
temperature  also  decreases,  but  a  little  more  rapidly  due  to  the  combination  of  the 
negative  heat  flux  and  decreasing  6^.  At  equilibrium,  6^,  Me,  and  are  all 
unchanging,  hence  and  are  also  unchanging. 


Unlike  Wangara,  the  ocean  supplies  a  constant  source  of  moisture.  This  moisture  is 
readily  transported  upward  in  the  PEL  when  convection  is  strong.  This  causes  the  PEL 
to  moisten  with  time  (see  Figure  5.b.6).  This  causes  both  and  to  increase.  As 
a  increases  it  begins  to  have  an  impact  on  (q^)^  which  causes  to  increase  more 

slowly  until  a  decreases  again.  Just  as  for  the  potential  temperature,  the  variables  that  the 
updraft  and  downdraft  mixing  ratio  depend  on  are  unchanging  at  equilibrium,  thus  (q^)^ 

and  (q^  do  not  change  either. 

The  updraft  and  downdraft  properties  at  level  E  are  shown  in  Figures  7.e.6  and7.e.7. 
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Ocean  Experiment  Up/Down  Potential  Temperatures  at  Level  B 


Figure  7.e.6:  (d^)^  and  {0j)^  for  Ocean  Experiment. 


Ocean  Experiment  Up/Down  Mixing  Ratios  at  Level  B 
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Figure  7.e.7:  and  (qj)g  for  Ocean  Experiment. 

The  potential  temperature  increases  with  height  in  the  ocean  experiment.  The  updraft  and 
downdraft  potential  temperature  at  level  B  depend  on  changes  in  dg  assuming  the  flux 
contribution  is  small  compared  to  these  changes.  The  PBL  depth  increases  rapidly  during 
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the  early  portion  of  the  simulation.  This  causes  relatively  large  changes  in  6^  compared 
to  the  flux  contribution.  Also,  the  flux  contribution  is  small  initially  because  it  contains 
Me  in  the  denominator  which  is  large  for  about  the  first  1000  minutes.  Therefore,  the 
properties  are  largely  controlled  by  changes  in  the  PBL  depth.  Both  and  {9j)g 
increase  when  Ap^  increases,  and  they  decrease  when  decreases.  At  equilibrium, 
Ap^  is  unchanging  so  and  {dj)g  are  unchanging  as  well. 

The  mixing  ratio  decreases  with  height  in  the  ocean  experiment,  but  the  ocean 
moistens  the  PBL  through  convection.  The  moistening  dominates  over  drying  that  occurs 
due  to  ascent.  Therefore,  increases  which  cause  (q^)^  and  {qj)g  to  increase  until  the 

mixed  layer  mixing  ratio  reaches  equilibrium.  At  this  point  q^  no  longer  changes. 

7.f.  Ocean  Experiment  PBL  Interior  Results 

The  ocean  experiment  interior  results  were  done  in  the  same  manner  as  Wangara 
using  a  f  of  1  second.  These  results  are  also  f  times  their  values  at  f  =  1  second  for 
other  settings  of  f .  The  dissipation  time  scale  is  shown  in  Figure  7.f.l. 
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Figure  7.f.l:  for  Ocean  Experiment. 
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The  two  peaks  correspond  to  a=l/2  (Tj^  For  Wangara  the  minimum  in 

occurred  when  a«l.  In  this  case  the  minimum  occurs  for  a  =  0.95.  During  the 
convective  period  is  about  2  orders  of  magnitude  longer  (not  considering  the  peaks) 

than  for  Wangara.  This  would  indicate  that  dissipation  was  more  effective  for  the 
Wangara  simulation  due  to  the  intense  convection. 

The  next  set  of  figures  contain  the  dissipation  rates  for  0  and  q  at  levels  S  and  B. 


Ocean  Experiment  Potential  Temperature  Dissipation  Rates 


Time  (minutes) 

Figure  7.f.2:  and  (eg)^  for  Ocean  Experiment. 
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Ocean  Experiment  Mixing  Ratio  Dissipation  Rates 


Time  (minutes) 

Figure  7.f.3:  and  for  Ocean  Experiment. 

The  initial  peaks  in  the  potential  temperature  dissipation  rates  are  predominantly  due  to 
the  surface  heat  flux.  The  rates  go  to  zero  when  a=l/2.  The  second  peak  in  is 

caused  by  a  large  negative  surface  heat  flux  and  a  minimum  in  the  TKE.  The  minimum 
in  (£g)g  that  occurs  at  the  same  time  is  caused  by  a  minimum  in  the  entrainment  rate.  At 

equilibrium,  dissipation  is  dominated  by  entrainment.  The  inversion  maintains  a 
temperature  gradient  at  the  top  of  the  PEL  which  creates  a  downward  flux.  The  small 
negative  heat  flux  at  the  surface  results  in  a  smaller  value  of  . 

The  large  initial  surface  moisture  flux  creates  the  first  peak  in  The  second 

peak  is  due  to  a  minimum  in  the  TKE  and  a  relatively  large  surface  flux.  The  minimum 
in  at  the  same  time  is  caused  by  a  minimum  in  E.  In  equilibrium,  the  dissipation 

rates  are  equal  because  the  surface  and  PEL  top  moisture  fluxes  are  equal.  The  moisture 
gradient  at  the  PEL  top  is  greater  than  at  the  surface,  but  V>E. 

The  last  diagnostics  for  the  ocean  experiment  are  the  potential  temperature  and 
mixing  ratio  gradients.  Figure  7.f.4  shows  the  potential  gradients  for  f  equal  to  10,  100, 
and  1000  seconds.  Like  Wangara,  the  I  second  gradients  were  too  large  and  are  not 
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shown.  The  mixing  ratio  gradients  for  100  and  1000  seconds  are  in  Figure  7.f.5.  The  10 
second  gradient  was  also  too  large. 


Ocean  Experiment  Potential  Temperature  Gradient 


Time  (minutes) 

dd 

Figure  7.f.4:  —  with  f  =  10, 100,  and  1000  Seconds  for  Ocean  Experiment. 
az 


Ocean  Experiment  Mixing  Ratio  Gradient 


Time  (minutes) 

Figure  7.f.5:  ^  with  f  =  10, 100,  and  1000  Seconds  for  Ocean  Experiment. 
az 
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All  the  gradients  show  the  potential  temperature  increasing  with  height  and  the  mixing 
ratio  decreasing  with  height,  except  at  the  very  beginning  of  the  simulation.  The  10 
second  gradients  appear  to  be  too  large  as  was  found  for  Wangara.  A  f  between  100  and 
1000  seconds  seems  most  suitable  for  this  type  of  simulation  as  well. 
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8.  Description  of  Two-Layer  Model 

The  two-layer  model  uses  the  same  set  of  equations  and  the  same  parameterizations  as 
the  one-layer  model  except  for  the  mixed  layer  potential  temperature  and  mixing  ratio 
equations.  Infinitesimal  ventilation  (surface)  and  entrainment  layers  are  maintained  with 
the  top  of  the  ventilation  layer  still  at  level  S  and  the  bottom  of  the  entrainment  layer  still 
at  level  B.  The  mixed  layer,  however,  is  divided  into  2  layers.  Level  1  is  within  the  top 
layer  and  level  2  is  in  the  bottom  layer.  The  layers  are  divided  at  level  I  (interior).  Figure 
8.1  is  a  diagram  of  the  two-layer  model. 


Entrainment  Layer 


Surface  Layer 


Figure  8.1:  Illustration  of  2-Layer  Model. 

Level  2  was  set  1/4  of  the  way  up  in  the  mixed  layer,  level  I  in  the  center,  and  level  1 
3/4  of  the  way  up.  The  levels  are  evenly  spaced  for  mathematical  ease.  The  equations 
used  to  predict  the  mean  potential  temperature  and  mixing  ratio  at  levels  1  and  2  do  not 
require  the  levels  to  be  equally  spaced.  These  equations  are  developed  in  the  next 
section.  Once  the  mean  potential  temperature  and  mixing  ratio  are  initialized  or  predicted 
at  levels  1  and  2,  the  mixed  layer  values  are  determined  using 


107 


(8.1) 


and 


e 

m  2 


Qn, 


(I1+Q2 

2 


(8.2) 


8.a.  Two-Layer  Potential  Temperature  and  Mixing  Ratio  Equations 

From  Randall  (personal  communication,  1993),  the  two-layer  equations  for  the 
potential  temperature  and  mixing  ratio  are 


^[{P,  -  Pb)^,]  =  -V  •[v;(p,  - Pb)^]  +  Ap„k,0,  + 

8[{Fs),+EeB,] 


(8.a.l) 


at  level  1  and 

|■[(Fs  -  Pi%]  =  -V  •[v2(p5  -  P,)^,]  -  + 

^  (8.a.2) 

at  level  2,  and 

i[{p,  -  Pb%]  =  -V  •  [v,{p,  -  + 

^  (8.a.3) 

at  level  1  and 

^[{Ps-  Pi%]  =  -'^  *[v2{Ps-  P,)^]- + 

^  (8.a.4) 

at  level  2.  In  these  equations,  is  the  vertical  velocity  at  level  I  as  seen  following  the 

coordinate  where  =  — — — . 

^Pm 
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Tlie  4  coordinate  is  similar  to  the  modified  (j-coordinate  used  by  Suarez  et  al. 
(1983).  The  coordinate  system  is  designed  so  that  the  earth's  surface  and  PBL  top  are 
coordinate  surfaces.  At  the  earth's  surface,  ^  =  0,  and  at  the  PBL  top,  ^  =  1.  For 

p.  <  p<  Pg,  ^  The  vertical  velocity,  ^ ,  measures  how  fast  a  ^  surface  moves 

4p„ 


as  the  PBL  depth  changes,  and  can  be  given  by  pressure 


final  pressure 


/initii 


initial  pressure 


At 


At  the  earth's  surface,  ^  is  always  0  no  matter  how  much  the  PBL  depth  changes,  so  is 


0  here.  At  the  PBL  top,  ^  is  always  1,  and  depends  on  how  much  the  PBL  depth  has 
changed. 


A  more  useful  formula  for  ^  is  obtained  by  adding  the  mass  conservation  equations 

3  — 

for  layers  1  and  2,  —{p,  -pg)  =  -V  •[v;(p,  -  pg)]  +  ApJ,  +  gE  and 

d  — 

^  (Ps  -  P/ )  =  (P5  -Pi)]-  4p.^/  .  together  to  get 

+  =  (8.a.5) 


and  then  using  the  conservation  of  mass  for  layer  2,  (8.a.5),  and  = 


Ps-Pi 

4p« 


to  obtain 


^Ji  =  ^/ V  •  -  V2  )4p„  ]  -  ^igE  ■ 


The  vertical  velocity  is  simply 


(8.a.6) 


Equations  (8.a.l)-(8.a.4)  can  be  written  in  advective  form  by  using  the  conservation 
of  mass  equations  for  the  two  layers.  Then,  assuming  horizontal  homogeneity,  except  for 
the  mean  divergence,  gives 

3ft  —  —  —  r  _  _  1 

(ft  =  -»-)]• 
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{Ps--P,)^  =  {^1  -^2)  +  )5  ~  )/ ] ’ 

(8.a.8) 

(P/  ~Pb)^  =  {d,  -q,)  +  -?;)]’ 

(8.a.9) 

(ps  -p/)^=4p.i,(^/ 

(8.a.l0) 

and  t,  =  • 

'  4p. 

(8.a.ll) 

Equations  (8.a.7)-(8.a.l0)  can  be  solved  if  the  flux  of  the  potential  temperature  and 
mixing  ratio  at  level  I  are  known.  The  equations  for  this  flux  are  developed  in  the  next 
section.  The  mean  value  of  6  and  q  at  level  I  is  just  the  mixed  layer  value  of  these 
variables  (0„=0i,  q^-  ^,).  Equations  (8.a.7)-(8.a.l0)  can  then  be  rewritten  in  terms  of 
6,  ,  q, ,  6, ,  and  ^3  using  (8.1),  (8.2),  and  the  interpolation  relations  for  and 

QsandB  ‘ 


The  mean  potential  temperature  increases  linearly  with  height,  and  the  mean  mixing 
ratio  decreases  linearly  with  height.  The  mean  value  of  these  variables  at  any  pressure  is 
given  by 

0(p)  =  a  +  6p  (8.a.l2) 

and 

q{p)  =  c  +  dp  (8.a.l3) 

Equations  (8.a.l2)  and  (8.a.l3)  are  just  equations  for  lines  where  a  and  c  are  intercepts, 
and  b  and  d  are  slopes  of  the  lines.  The  slopes  are  given  by  the  difference  of  the  mean 
values  at  levels  1  and  2  divided  by  the  difference  in  pressure  between  levels  1  and  2.  The 
intercepts  are  then 


0(p  =  l)  =  ».=a  +  ^l-^p, 

P,-P2 


or 


a  =  6^ 


e,-e, 

P1-P2 


Pn 


(8.a.l4) 
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and 


q(p  =  I)  =  q^=c  +  ^J—^p, 
P1-P2 


or 


-^P,. 

P1-P2 


(8.a.l5) 


Inserting  (8.a.l4)  into  (8.a.l2)  and  (8.a.l5)  into  (8.a.l3)  gives  the  mean  quantities  at  any 
pressure, 

P)  =  ^m  +  {P-Pi) (8.a.l6) 
Pi  -P2 

and 

^(P)  =  Qn,+{p-Pi)--^^-  (8.a.l7) 

Pi  -  P2 


The  interpolation  relations  for  0  and  q  at  levels  S  and  B  are  obtained  by  using  p=ps  and 
P=Pb  in  (8.a.l6)  and  (8. a.  17). 

Then,  the  finite  difference  forms  of  (8.a.7)-(8.a.l0)  for  Wangara  using  a  backward 
(implicit)  scheme  are 


er  - 


gE/M(  Ps-p, 
^Pn,  [Pi-Pb 


^  .  gEAt^ 


0r  + 


4p„ 


Ps-Pi 

\P1-P2) 


0^  - 


2^{Pi-P2) 


4?, 


m  \ 


Ps-Pi 

PI-P2J 


+ 


{Ps-p,)o{I-o)^t^„,i  ^  gEAtf  Ps-P,^ 
2i{Pi-P2)  '  ^n.[  Pi-Pb 


K+^{F.)s 


^n. 


{Pi-Pb) 


(8.a.l8) 
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nn*l  gEAt 


Ap„ 


■e; 


Ap„ 


P,-Pb 

KPi-P:J 


]dr‘  I 


2i{Pj-P:) 


j^At  ^gEAt^ 

k  i  * 


{Pi-Pb) 


Ap^ 


P,~Pb 

\Pi~Pzj 


/«  +  / 


{P,-PBW-(^)^t-,,  _  gAt  ,^,_gEAt^ 

2i{p,-P:)  ■  ~{Ps-P.r'^^  Ap„^~ 


gAt 


'  P,-Pb^ 
yPs-Pi) 


{Eb\-^0^.. 


and 


Q, 


(  Ps-pY 

iPs-P,^ 

[pi-PbJ 

Ap^ 

[P,-Pzj 

,'—/»+/ 

■ 


{Ps-p,)<y(}-<r)At^„^,  _jE^  ,  ^sY^-n.,  ^ 
2i{p,-P:)  '  {Pi-Pb)^‘  2Ap^ 


gVAt 


f 


Ap„ 


Ps-P, 


VPj-P: 


i~n+/ 


gEAt 


z  J 


Ap„ 


\ 


Ps-P, 
KP,-Pzj 


I— n+/  ,  {Ps-P,)o{l-o)At 

— n+/ 

p-  ^  \ 


2i{p,-Pz) 


24p.  ■ 


;t  4.77'' 
■<7fi+  Qt  ’ 


gVAt 

(ps-P,) 

_  gEAt 

(  \ 
Ps  -  P, 

Ap„ 

[Pj-PzJ 

'  ■  Ap„, 

[p,-Pb) 

gVA(_ 

QB*+-T~(is- 


{p,-Pb) 


gEAt  ,  gEAt 

40."' 


4p. 


'  P,-Pb^ 
yP,-Pz) 


^7”*^  4-  — /!+/  , 


P,  -  Pfi 


2{Ps-P,)  ‘  {P,-Pz) 


gVAt 


Ap„ 


P,-Pb 
{  P,  -  Pz  J 


+  SE  At 


2Ai?„, 


\Ps-P, 


^n, 


P,-Pb 
P,  -  Pz  ) 


^r‘  - 


{p,-PB)(^il-Cf)At  — /!+/  ,  gVAt 

— n+/  gVAt 

— 

2i{p,-Pz)  '  2{ps~p,f'  {P,-Pzf' 


gVAt 


f  \ 

P,  -Pb 

i  P,~Pb] 

V  Ps  -  P,  ) 

Al^. 

1  P,  ~  Pz  J 

2AtK 


gEAt  _  gVAt 

■^fl+ 


gV.4/  _ 

+ 


{Ps-Pi) 


Ap„. 


Ap„ 


l  Ps  -  P,  j 


(8.a.l9) 


+ 

(8.a.20) 

+ 


(8.a.21) 
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The  surface  heat  flux  appears  explicitly  in  (8.a.l8)  and  (8.a.l9)  because  it  was  prescribed 
for  the  Wangara  simulations.  For  the  ocean  simulations  (8.a.20)  and  (8.a.21)  were  used 
with  6  in  place  of  q,  instead  of  (8.a.l8)  and  (8.a.l9),  for  the  potential  temperatures  at 
levels  1  and  2. 

Equations  (8.a.l8)  and  (8.a.l9)  are  two  equations  in  two  unknowns,  6j  and  0,. 
Equations  (8.a.20)  and  (8.a.21)  are  also  two  equations  in  two  unknowns,  q,  and  q^. 

These  sets  of  equations  are  solved  simultaneously  to  obtain  the  mean  values  at  levels  1 
and  2.  The  mixed  layer  values  are  finally  determined  using  (8.1)  and  (8.2). 

8.b.  Two-Layer  Model  Diagnostics 

Since  the  two-layer  model  predicts  the  mean  values  of  0  and  q  at  levels  1  and  2,  the 
gradients  of  these  variables  were  determined  by  using 


and 


II 

\  t 

(8.b.l) 

dq  _  qj-q. 
dz 

(8.b.2) 

instead  of  equation  (6.d.9). 

The  final  form  for  the  flux  is  obtained  by  truncating  ((6.d.l  1)  at  first  order  in 
I  dp,  and  substituting  this  into  (6.d.l0)  at  level  1.  This  gives 


'Ps-Pl^ 

\  J 


■('v),, 


Pi  Pb 


-I- 


HPi  -  Pb){Ps  -  PfW dy/ 
2  gr  dp 


(8.b.3) 


The  gradient  in  (8.b.3)  is  determined  using  (8.b.l)  or  (8.b.2)  and  the  hydrostatic  relation. 
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9.  Two-Layer  Model  Results 


Prognostic  results  for  the  two-layer  model  using  the  Wangara  data  are  presented  and 
compared  with  the  prognostic  results  from  the  one-layer  model.  Diagnostic  results  for  the 
gradients,  and  the  mean  values  of  0  and  q  at  levels  S,  2, 1, 1,  and  B  are  also  shown.  The 
diagnostic  results  were  obtained  using  a  f  of  1, 10, 100,  and  1000  seconds. 

9.a.  Two-Layer  Prognostic  Results 

Figures  9.a.l-9.a.5  show  the  prognostic  variables  using  the  two-layer  model. 


0  240  480  720  960  1200  1440 

Time  (minutes) 


Figure  9.a.l:  Two- Layer 
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Two-Layer  Potential  Temperature 


Time  (minutes) 
Figure  9.a.2:  Two-Layer 


Two-Layer  Mixing  Ratio 


Time  (minutes) 


Figure  9.a.3:  Two-Layer 


10 


Two-Layer  Horizontal  Velocity 


Tinie  (minutes) 
Figure  9.a.4:  Two-Layer  |v„|. 


Two-Layer  Turbulence  Kinetic  Energy 


Time  (minutes) 

Figure  9.a.5:  Two-Layer 

This  is  no  difference  between  these  figures  and  the  figures  for  the  one-layer  prognostic 
variables.  This  was  expected  for  Az^,  |v„|,  and  e„  which  are  predicted  the  same  way  in 

each  model.  Identical  values  for  and  indicate  that  the  two-layer  model  is 


116 


functioning  properly.  Prognostic  variables  for  the  two-layer  ocean  experiment  are  not 
shown,  but  they  were  also  the  same  as  the  one-layer  ocean  experiment  variables. 

9.b.  Two-Layer  Diagnostic  Results 

To  show  the  effects  of  varying  f ,  the  initial  values  of  the  mean  6  and  q  were  set  to 

initio  ~ 

(^2 

initial  ~  i^m)ininal  ~  ^  S  ^8  ' 

and 

i^2Ltial  =  i^Xtia,  +  ^  8  k8'‘-  (9.b.l) 

This  gave  a  sounding  where  the  potential  temperature  was  initially  increasing  with  height 
and  the  mixing  ratio  was  initially  decreasing  with  height.  The  gradients  of  d  and  q  were 
then  determined  with  these  initial  conditions  and  are  shown  in  Figures  9.b.l  and  9.b.2. 


Two-Layer  Potential  Temperature  Gradient 
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Figure  9.b.l  :  —  with  T  =  1, 10, 100,  and  1000  Seconds  for  Two-Layer  Model. 
az 
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Two-Layer  Mixing  Ratio  Gradient 


0  10  20  30  40  50  60 

Time  (minutes) 

Figure  9.b.2:  ^  with  f  =  1, 10, 100,  and  1000  Seconds  for  Two-Layer  Model. 
az 

The  gradients  for  f  of  1  and  10  seconds  were  very  small  all  the  time.  The  1000  second 
gradients  are  the  only  ones  that  were  not  near  zero  within  60  minutes.  The  100  second 
gradients  started  out  steep,  but  did  decrease  to  near  zero  by  60  minutes.  These  gradients 
were  created  by  the  artificial  initial  conditions  in  the  mean  values  of  8  and  q  at  levels  1 
and  2.  The  actual  gradients  for  all  ts  were  near  zero.  This  and  the  one-layer  gradients  do 
indicate,  however,  that  a  f  not  much  larger  than  100  seconds  should  be  used. 

Figures  9.b.3-9.b.6  depict  the  mean  potential  temperature  soundings  using  Q  at 
levels  S,  2.  1,  and  B,  for  the  different  values  of  f  . 
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Two-Layer  Mean  Potential  Temperatures  For  1  Second 


Figure  9.b.3:  Mean  Potential  Temperatures  Using  f  =  1  Second. 


Two-Layer  Mean  Potential  Temperatures  For  10  Seconds 


Figure  9.b.4:  Mean  Potential  Temperatures  Using  f  =  10  Seconds. 
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Two-Layer  Mean  Potential  Temperatures  For  100  Seconds 


Time  (minutes) 


Figure  9.b.5;  Mean  Potential  Temperatures  Using  f  =  100  Seconds. 


Two-Layer  Mean  Potential  Temperatures  For  1000  Seconds 


Time  (minutes) 

Figure  9.b.6:  Mean  Potential  Temperatures  Using  f  =  1000  Seconds. 

These  figures  indicate  how  long  it  took  for  the  mean  potential  temperatures  to  adjust  to 
the  mixed  layer  potential  temperature.  For  f  =  1  second,  it  only  took  3  minutes  for 
adjustment.  When  f  was  set  to  10  seconds,  the  adjustment  time  increased  to  about  15 
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minutes.  At  a  i  of  100  seconds,  the  time  to  adjust  had  jumped  to  a  little  over  an  hour. 
Also,  the  temperatures  diverged  for  short  periods  twice.  Finally,  when  f  was  set  to  1000 
seconds,  the  temperatures  did  not  adjust  until  t=700  minutes  (12  hours).  There  was 
considerable  divergence  in  the  temperatures  initially,  and  a  small  amount  of  divergence 
from  about  t=240  to  300  minutes. 
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10.  Summary  and  Conclusions 

A  single-layer  bulk  boundary  layer  model  was  presented  that  predicts  the  mixed  layer 
values  of  the  potential  temperature,  mixing  ratio,  and  u  and  v  momentum.  The  model  also 
predicts  the  depth  of  the  boundary  layer  in  terms  of  pressure  (Apm)  and  the  turbulence 
kinetic  energy  (TKE).  The  TKE  prediction  equation  was  formulated  using  a  second-order 
closure  that  relates  the  dissipation  velocity  to  the  TKE.  The  model  also  diagnostically 
determines  the  fractional  area  covered  by  rising  motion  (o)  and  the  entrainment  rate  (E). 

Positive  and  negative  entrainment  rate  parameterizations  were  developed,  and  the  one 
used  for  a  particular  time  step  was  based  on  the  sums  of  the  buoyancy  (B)  and  shear  (S) 
production  (with  and  without  E  included).  A  tunable  parameter  was  used  to  specify  a 
fraction  of  the  sums  to  check.  This  was  done  to  prevent  a  large  positive  E  from  suddenly 
becoming  negative.  A  value  of  0.9  for  this  parameter  was  found  to  produce  good  results. 

The  positive  entrainment  rate  was  parameterized  by  assuming  that  E  is  proportional  to 
the  square  root  of  the  TKE.  The  constants  in  the  parameterization  were  obtained  by 
assuming  a  balance  between  buoyant  production  and  dissipation,  and  using  large-eddy 
simulation  results  from  Deardorff  (1974).  This  parameterization  led  to  two  Richardson 
number  limits,  Ri»l  (strong  inversion)  and  Ri=0  (no  inversion). 

The  negative  entrainment  rate  was  parameterized  by  assuming  that  E  and  em  are  small 

compared  to  their  values  during  rapid  PBL  growth.  The  local  change  term  was  then 

neglected  in  the  Cm  equation  which  led  to  a  balance  between  the  entrainment  rate  and 
B  +  S  ^  D 

- .  A  tunable  parameter  was  then  introduced  to  partition  this  balance  equation  into 

a  weighted  contribution  of  the  local  change  of  em  and  the  production  of  Cm  due  to  E.  A 
value  of  0.9  was  used  for  the  simulations  and  produced  the  best  results. 

Two  simulations  were  run.  The  first  simulation  used  the  Wangara  Day  33  PBL  data. 
The  surface  heat  flux  was  prescribed  using  a  sine  approximation.  The  ventilation  (surface) 
mass  flux  was  parameterized  using  the  formulation  from  Louis  (1979)  and  was  used  for  the 
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surface  momentum  and  moisture  fluxes.  The  land  simulation  was  initialized  using  the 
Wangara  data. 

The  diurnal  trend  of  the  mixed  layer  depth,  except  for  night  values,  was  accurately 
depicted  by  the  model.  The  model  captured  the  slow  growth  early  in  the  morning  when 
there  was  a  strong  inversion,  rapid  growth  during  mid-morning  when  the  inversion  broke, 
slow  growth  during  the  afternoon  under  a  quasi  steady-state  PBL  topped  by  an  inversion, 
and  rapid  decay  after  loss  of  surface  heating  at  sunset.  The  nocturnal  PBL  did  not  grow 
slowly  as  expected.  There  appears  to  be  a  problem  with  the  negative  entrainment 
parameterization  at  night.  The  shear  production  at  night  due  to  the  nocturnal  jet  should  be 
sufficient  to  allow  the  PBL  to  grow  even  with  negative  buoyancy  production. 

Diagnostic  variables  to  study  the  characteristics  of  a  clear  convective  boundary  layer 
(CBL)  were  developed  using  the  concept  of  the  convective  mass  flux  model.  Equations 
were  presented  for  the  plume-scale  variance  transport  of  a  scalar,  y,  and  updraft  and 
downdraft  properties  of  y.  Then  the  convective  mass  flux  was  matched  with  the 
ventilation  and  entrainment  layer  fluxes.  This  was  accomplished  by  assuming  these  layers 
were  infinitesimal,  and  the  small-eddy  fluxes  at  levels  S  and  B  were  negligible  compared 
with  the  convective  circulations.  Use  of  the  TKE  then  allowed  the  convective  mass  flux 
and  the  fractional  area  covered  by  rising  motion  to  be  determined  using  model  variables. 

The  features  of  the  CBL  were  well  illustrated  by  the  model  diagnostic  results.  The 
model  showed  the  dominance  of  buoyancy  production  over  shear  production  in  a  CBL. 
This  was  shown  by  a  plot  of  the  buoyancy  production  versus  the  shear  production,  and  by 
a  plot  of  the  plume-scale  variance  transport  of  9  at  levels  S  and  B.  The  entrainment  rate 
was  also  shown  to  be  an  important  mechanism,  especially  during  rapid  growth  when  E 
became  large.  The  intense  convection  typical  of  a  CBL  was  indicated  by  o«l.  The 
convective  mass  flux  was  a  minimum  at  this  time,  contrary  to  what  one  would  expect. 
However,  during  vigorous  convection  when  E»V,  (y«l,  and  Me  is  small  because 

The  updraft  and  downdraft  properties  further  highlighted  the  CBL  characteristics.  The 
updrafts  at  level  S  were  warmer  and  wetter  than  the  dow'ndrafts.  Here,  the  convection  was 
seen  in  terms  of  the  surface  heating  rate  and  the  surface  heat  and  moisture  transport  rates. 
The  dominance  of  one  of  these  over  the  other  was  important  in  determining  the  behavior  of 
the  updraft  and  downdraft  properties  at  the  surface. 
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The  inversion  at  the  top  of  the  PBL  was  the  controlling  factor  for  the  updraft  and 
downdraft  properties  at  level  B.  When  there  was  an  inversion,  the  updraft  air  was  cooler 
and  wetter  than  the  downdraft  air.  When  the  inversion  was  absent,  the  updraft  air  was 
warmer  and  much  wetter  than  the  downdraft  air.  This  was  caused  by  the  strong  convection 
that  rapidly  transported  heat  and  moisture  upwards. 

Diagnostics  for  the  PBL  interior  were  developed  to  gain  further  insight  into  the  CBL. 

A  balance  equation  was  presented  for  the  variance  of  \j/.  Each  term  in  this  equation  was 
modeled  to  obtain  equations  for  the  variance  and  dissipation  rate.  A  dissipation  time  scale 
in  terms  of  the  model  parameter  f  was  introduced.  The  balance  equation  was  then  solved 
to  get  a  relation  for  the  gradient  of  ij/. 

The  dissipation  time  scale  was  found  to  be  the  shortest  during  the  period  when  the 

surface  heating  was  the  strongest,  corresponding  to  the  high  efficiency  of  the  small-eddies. 
As  expected  at  this  time,  the  dissipation  rates,  £g  and  e^,  were  at  their  largest  values.  The 

dissipation  rate  of  0  at  the  surface  dominated  over  the  dissipation  rate  at  level  B.  Again, 
this  was  due  to  the  strong  surface  heating  present.  For  q,  the  opposite  was  true.  The  large 
value  of  E  caused  the  moisture  flux  at  the  PBL  top  to  be  much  greater  than  at  the  surface, 
especially  since  mixing  had  reduced  the  surface  to  mixed  layer  moisture  gradient. 

The  gradients  of  0  and  q  were  determined  using  a  f  of  10,  100,  and  1000  seconds. 
The  gradient  results  were  matched  to  the  Wangara  data  to  determine  the  best  value  for  f.  A 
value  between  100  and  1000  seconds  seemed  most  reasonable  based  on  the  data.  The  1000 
second  gradients  showed  the  expected  increase  in  potential  temperature  with  height  and 
decrease  of  moisture  with  height,  typical  of  a  convective  boundary  layer. 

A  surface  transfer  coefficient  was  developed  using  the  TKE,  and  was  determined  to  be 
about  an  order  of  magnitude  larger  than  the  transfer  coefficient  normally  found  in  the  bulk 
aerodynamic  formula  for  V.  This  was  expected  because  the  surface  velocity  was  about  10 
times  the  square  root  of  the  TKE.  Using  this  transfer  coefficient  over  the  conventional  one 
has  the  advantage  that  V  exists  if  there  is  turbulence,  even  if  the  surface  wind  is  zero.  This 
may  occur  in  a  heated  boundary  layer  where  turbulence  is  generated  only  by  buoyancy 
when  the  surface  wind  is  calm. 
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The  period  when  the  inversion  vanished  was  clearly  indicated  by  the  Richardson 
number.  The  limit  for  Ri=0  was  about  1  during  this  time  as  expected.  When  the  inversion 
was  strong  in  the  afternoon,  the  limit  for  Ri»l  approached  0.2,  but  was  too  small.  The 
assumption  of  the  balance  between  buoyant  production  and  dissipation  that  led  to  the 
relation  for  the  limit  when  Ri»l  may  be  slightly  inaccurate. 

A  one-layer  simulation  using  simple  ocean  data  was  then  mn  to  obtain  steady-state 
solutions.  Fixed  surface  and  top  mixing  ratios,  sea  surface  temperature,  surface  winds, 
and  geostrophic  winds  were  used.  The  temperature  and  winds  at  the  top  of  the  PBL  were 
determined  by  constant  lapse  rates.  The  surface  fluxes  of  heat,  moisture,  and  momentum 
were  determined  using  Louis  (1979)  ventilation  mass  flux  formulation.  A  divergence  of  4 
10-6  5-1  ^as  used  to  balance  E  in  the  Apm  prediction  equation. 

The  prognostic  variables  converged  to  their  equilibrium  values  by  100  hours.  The 
steady-state  form  of  the  prognostic  equation  for  0  was  derived.  This  equation  was  used  to 
compare  the  value  of  0  with  the  model  predicted  value.  The  value  from  the  equation  was 
only  0.04  K  different  from  that  predicted  by  the  model. 

In  a  study  done  by  Schubert  et  al.  (1979),  they  found  the  adjustment  time  for  the  PBL 

depth  was  considerably  longer  than  for  the  other  prognostic  variables  when  the  ratio  of 
C  V 

- - - was  about  4  or  5.  This  would  indicate  that  surface  transfer  dominates  over 

Dzg  -hdZgldt 

mixing  at  the  PBL  top.  The  value  of  the  ratio  obtained  in  the  present  ocean  simulation  was 
not  constant,  but  was  never  less  than  4.5.  However,  the  adjustment  time  of  the  PBL  depth 
was  the  same  compared  to  the  other  prognostic  variables. 

This  discrepancy  can  be  explained  by  tlie  different  enuninment  parameterization  used  in 
Schubert's  study  and  the  present  model.  In  Schubert's  study  a  constant  small  value  of  E 
was  used,  while  in  the  present  study  E  varied  and  became  large  during  the  early  portion  of 
the  simulation.  The  large  value  of  E  allowed  the  PBL  depth  to  adjust  as  fast  as  the  other 
prognostic  variables.  This  was  shown  by  determining  how  long  adjustment  would  have 
taken,  had  E  been  small  and  constant  during  the  entire  simulation.  This  adjustment  time 
was  about  77  days,  which  corresponds  to  an  adjustment  time  for  the  PBL  depth  of  about 
20  times  as  long  as  the  adjustment  for  the  other  variables.  This  agrees  with  the  results 
obtained  in  Schubert's  study. 
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A  two-layer  model  that  predicts  the  mean  values  of  0  and  q  at  two  levels  in  the  PBL 
was  then  developed  to  address  the  problem  of  the  large  gradients  obtained  by  the  one-layer 
model.  The  model  was  developed  by  equally  spacing  the  levels  for  mathematical 
simplicity,  even  though  the  2-level  equations  do  not  require  these  constraints.  This  model 
retains  all  the  parameterizations  used  in  the  one-layer  model.  The  only  differences  are  the 
determination  of  the  mixed  layer  values  of  0  and  q,  and  the  gradients  of  6  and  q . 

The  two-layer  model  produced  the  same  results  for  the  prognostic  variables  using  the 
Wangara  data  as  the  one-layer  model.  This  verified  that  the  model  worked  correctly.  The 
gradients  of  6  and  q  were  near  zero  for  the  entire  simulation  which  differed  considerably 
from  the  one-layer  model  gradients.  Identical  results  were  also  obtained  for  the  ocean 
experiment. 

TTie  initial  values  of  the  mean  values  of  0  and  q  at  levels  1  and  2  were  perturbed  to 
study  the  effects  of  changing  f .  The  gradients  were  found  to  be  larger  at  a  given  time  step 
as  f  was  increased.  The  gradients  for  all  values  of  f  except  1 000  seconds  approached 
zero  within  60  minutes.  Also,  the  mean  values  of  the  potential  temperatures  at  levels  S,  2, 
1,  and  B  converged  to  the  mixed  layer  potential  temperature  within  60  minutes  for  all 
values  of  f  except  1000  seconds.  This  result,  along  with  the  gradients  from  the  one-layer 
model,  indicate  that  a  f  near  100  seconds  is  the  best  choice. 

Following  is  a  summary  of  items  that  were  presented  for  the  first  time  in  this  thesis: 

(1)  A  positive  entrainment  rate  parameterization  that  assumed  a  balance  between 
buoyancy  production  and  dissipation  of  turbulence  kinetic  energy. 

(2)  A  negative  entrainment  rate  parameterization  tliat  allowed  the  PBL  depth  to  decrease 
late  in  the  day  when  buoyancy  production  was  no  longer  sufficient  to  maintain  the 
turbulence. 

(3)  A  fully  implicit  finite  difference  equation  for  the  TKE  (when  the  entrainment  rate  is 
positive)  solved  as  a  cubic  equation.  The  square  of  the  solution  that  is  always  real  was 
assigned  to  the  TKE. 

(4)  Results  for  both  the  Wangara  and  Ocean  studies  showing  the  fractional  area 
covered  by  rising  motion,  convective  mass  flux,  updraft  and  downdraft  properties  of  6 
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and  q  at  the  surface  and  PBL  top,  dissipation  rates  of  0  and  q  at  the  surface  and  PBL  top. 
dissipation  time  scale,  and  gradients  of  d  and  q . 

(5)  Results  and  comparison  for  the  Wangara  study  of  two  surface  bulk  transfer 
coefficients,  one  dependent  on  the  surface  velocity  and  the  other  on  the  turbulence  kinetic 
energy. 

(6)  A  two-layer  model  which  predicted  6  and  q  at  two  levels. 

(7)  Equations  that  determined  the  upward  turbulent  fluxes  of  0  and  q  in  the  interior  of 
the  PBL.  These  equations  were  used  to  obtain  6  and  q  in  the  two-layer  model. 

The  one  and  two-layer  models  presented  provide  an  accurate  representation  of  the  clear 
CBL.  The  turbulence  characteristics  are  depicted  by  the  prognostic  turbulence  kinetic 
energy  equation.  However,  the  PBL  typically  contains  clouds.  Future  work  should 
include  adding  cloud  effects  to  these  models.  TTiis  can  be  approached  in  two  steps.  First, 
a  simplified  dry  cloud  layer  should  be  added  which  would  have  the  effect  of  radiatively 
cooling  the  air  above  the  cloud.  This  is  a  relatively  simple  step.  Next,  as  a  more  complex 
procedure,  moist  processes  should  be  included.  Lilly  (1968)  provides  a  means  for 
accomplishing  these  steps. 

Additional  work  should  also  be  done  to  obtain  a  better  representation  for  the  nocturnal 
PBL.  The  positive  and  negative  entrainment  relations  would  have  to  be  modified.  The 
addition  of  mc.'e  complicated  radiative  processes  besides  a  simple  radiative  cooling  term, 
and  a  parameterization  that  takes  into  account  the  nocturnal  jet,  may  allow  the  PBL  to  grow 
at  night. 

The  convective  mass  flux  and  the  fractional  area  covered  by  rising  motion  were 
assumed  to  be  constant  with  height.  However,  large-eddy  simulations  indicate  that  these 
variables  are  not  constant  with  height.  Height  dependent  equations  for  these  variables 
should  be  developed.  Randall  et.  al.  (1992)  provides  a  possible  approach  to  accomplish 
this. 


The  two-layer  model  should  also  be  further  developed  with  the  above  suggestions.  In 
addition,  the  momentum  should  be  calculated  at  levels  1  and  2.  Then  all  the  prognostic 
variables  would  be  determined  at  the  same  resolution.  Next,  cloud  effects  should  be  added 
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to  allow  the  model  to  operate  in  a  wide  variety  of  meteorological  conditions.  The  model 
should  then  be  modified  to  make  predictions  at  multiple  levels.  Finally,  the  model  should 
be  incorporated  into  the  CSU  GCM. 
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